| Previous changeset 14:2efa46ce2c4c (2017-10-18) Next changeset 16:1710b0e874f1 (2017-10-21) |
|
Commit message:
fastqc_report v2.0.0 |
|
modified:
fastqc_report.Rmd fastqc_report.xml fastqc_report_render.R |
|
removed:
fastqc_report_ori.Rmd fastqc_report_render_ori.R |
| b |
| diff -r 2efa46ce2c4c -r d1d20f341632 fastqc_report.Rmd --- a/fastqc_report.Rmd Wed Oct 18 22:06:39 2017 -0400 +++ b/fastqc_report.Rmd Thu Oct 19 00:11:14 2017 -0400 |
| [ |
| @@ -1,5 +1,5 @@ --- -title: 'HTML report title' +title: 'Short reads evaluation with [FastQC](https://www.bioinformatics.babraham.ac.uk/projects/fastqc/)' output: html_document: number_sections: true @@ -35,18 +35,47 @@ done ``` -* Create links to original HTML reports +## Evaluation results ```{r 'html report links'} -html_report_list = list() -html_files = list.files('REPORT_DIR', pattern = '.*html') -for (i in html_files) { - html_report_list[[i]] = tags$li(tags$a(href=i, i)) -} -tags$ul(html_report_list) +html_file = list.files('REPORT_DIR', pattern = '.*html') +tags$ul(tags$a(href=html_file, paste0('HTML report', opt$name))) +``` + + +```{r 'extract fastqc_data.txt and summary.txt'} +# list all zip files +zip_file = list.files(path = 'REPORT_DIR', pattern = '.zip') +unzip(paste0('REPORT_DIR/', zip_file), exdir = 'REPORT_DIR') + +unzip_directory = paste0(tail(strsplit(opt$reads, '/')[[1]], 1), '_fastqc/') +fastqc_data_txt_path = paste0('REPORT_DIR/', unzip_directory, 'fastqc_data.txt') +summary_txt_path = paste0('REPORT_DIR/', unzip_directory, 'summary.txt') ``` -# Fastqc output summary + +```{r 'summary.txt'} +tags$ul(tags$a(href=paste0(unzip_directory, 'summary.txt'), 'summary.txt')) +``` + + +```{r 'fastqc_data.txt'} +tags$ul(tags$a(href=paste0(unzip_directory, 'fastqc_data.txt'), 'fastqc_data.txt')) +``` + + +# Fastqc output visualization + +## Overview + +```{r} +# read.table(fastqc_data_txt_path) +summary_txt = read.csv(summary_txt_path, header = FALSE, sep = '\t')[, 2:1] +names(summary_txt) = c('MODULE', 'PASS/FAIL') +knitr::kable(summary_txt) +``` + +## Summary by module {.tabset} * Define a function to extract outputs for each module from fastqc output @@ -62,7 +91,21 @@ } ``` -## +### Per base sequence quality + +```{r} +pbsq = extract_data_module(fastqc_data_txt_path, 'Per base sequence quality') +knitr::kable(pbsq) +``` + +### Per tile sequence quality + +```{r} +ptsq = extract_data_module(fastqc_data_txt_path, 'Per tile sequence quality') +knitr::kable(ptsq) +``` + + # Session Info |
| b |
| diff -r 2efa46ce2c4c -r d1d20f341632 fastqc_report.xml --- a/fastqc_report.xml Wed Oct 18 22:06:39 2017 -0400 +++ b/fastqc_report.xml Thu Oct 19 00:11:14 2017 -0400 |
| b |
| @@ -30,6 +30,7 @@ Rscript '${__tool_directory__}/fastqc_report_render.R' -e $echo -r $reads + -n $reads.name -o $report -d $report.files_path |
| b |
| diff -r 2efa46ce2c4c -r d1d20f341632 fastqc_report_ori.Rmd --- a/fastqc_report_ori.Rmd Wed Oct 18 22:06:39 2017 -0400 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 |
| [ |
| b'@@ -1,381 +0,0 @@\n----\n-title: "Fastqc report: short reads quality evaluation"\n-author: "Ming Chen"\n-output: html_document\n----\n-\n-```{r setup, include=FALSE}\n-knitr::opts_chunk$set(echo=ECHO, warning=FALSE, message=FALSE)\n-library(plyr)\n-library(stringr)\n-library(dplyr)\n-library(highcharter)\n-library(DT)\n-library(reshape2)\n-library(plotly)\n-library(formattable)\n-library(htmltools)\n-```\n-\n-\n-```{bash \'create output directory\', echo=FALSE}\n-# create extra files directory. very important!\n-mkdir REPORT_OUTPUT_DIR\n-```\n-\n-# Fastqc analysis\n-```{bash \'copy data to working directory\', echo=FALSE}\n-# Copy uploaded data to the working directory\n-for f in $(echo READS | sed "s/,/ /g")\n-do\n- cp $f ./\n-done\n-```\n-\n-\n-```{bash \'run fastqc\', echo=FALSE}\n-for r in $(ls *.dat)\n-do\n- fastqc -o REPORT_OUTPUT_DIR $r > /dev/null 2>&1\n-done\n-```\n-\n-## Fastqc html reports\n-\n-Below are links to ***Fastqc*** original html reports.\n-```{r \'html report links\'}\n-html_report_list = list()\n-html_files = list.files(\'REPORT_OUTPUT_DIR\', pattern = \'.*html\')\n-for (i in html_files) {\n- html_report_list[[i]] = tags$li(tags$a(href=i, i))\n-}\n-tags$ul(html_report_list)\n-```\n-\n-\n-## Parsing fastqc data\n-\n-```{bash echo=FALSE}\n-##==== copy fastqc generated zip files from report output directory to job work directory ==\n-cp -r REPORT_OUTPUT_DIR/*zip ./\n-\n-# create a file to store data file paths\n-echo "sample_id,file_path" > PWF_file_paths.txt # Pass, Warning, Fail\n-echo "sample_id,file_path" > PBQS_file_paths.txt # Per Base Quality Score\n-echo "sample_id,file_path" > PSQS_file_paths.txt # Per Sequence Quality Score\n-echo "sample_id,file_path" > PSGC_file_paths.txt # Per Sequence GC Content\n-echo "sample_id,file_path" > PBSC_file_paths.txt # Per Base Sequence Content\n-echo "sample_id,file_path" > PBNC_file_paths.txt # Per Base N Content\n-echo "sample_id,file_path" > SDL_file_paths.txt # Sequence Duplication Level\n-echo "sample_id,file_path" > SLD_file_paths.txt # Sequence Length Distribution\n-echo "sample_id,file_path" > KMC_file_paths.txt # Kmer Content\n-\n-for i in $(ls *.zip)\n-do\n- BASE=$(echo $i | sed \'s/\\(.*\\)\\.zip/\\1/g\')\n- echo $BASE\n- unzip ${BASE}.zip > /dev/null 2>&1\n- \n- ##====== pass,warning,fail (WSF) =============\n- awk \'/^>>/ {print}\' "$BASE"/fastqc_data.txt | grep -v \'END_MODULE\' | sed \'s/>>//\' > "$BASE"-PWF.txt\n- echo "${BASE},${BASE}-PWF.txt" >> PWF_file_paths.txt\n-\n- ##====== per base quality scores (PBQS) ======\n- awk \'/^>>Per base sequence quality/ {flag=1; next} /END_MODULE/ {flag=0} flag\' "$BASE"/fastqc_data.txt >"$BASE"-PBQS.txt\n- echo "${BASE},${BASE}-PBQS.txt" >> PBQS_file_paths.txt\n-\n- ##====== per sequence quality scores (PSQS)\n- awk \'/^>>Per sequence quality scores/ {flag=1; next} /END_MODULE/ {flag=0} flag\' "$BASE"/fastqc_data.txt >"$BASE"-PSQS.txt\n- echo "${BASE},${BASE}-PSQS.txt" >> PSQS_file_paths.txt\n-\n- ##====== Per sequence GC content (PSGC)\n- awk \'/^>>Per sequence GC content/ {flag=1; next} /END_MODULE/ {flag=0} flag\' "$BASE"/fastqc_data.txt >"$BASE"-PSGC.txt\n- echo "${BASE},${BASE}-PSGC.txt" >> PSGC_file_paths.txt\n- \n- ##====== Per Base Sequence Content (PBSC)\n- awk \'/^>>Per base sequence content/ {flag=1; next} /END_MODULE/ {flag=0} flag\' "$BASE"/fastqc_data.txt >"$BASE"-PBSC.txt\n- echo "${BASE},${BASE}-PBSC.txt" >> PBSC_file_paths.txt\n- \n- ##====== Per Base N Content (PBNC)\n- awk \'/^>>Per base N content/ {flag=1; next} /END_MODULE/ {flag=0} flag\' "$BASE"/fastqc_data.txt >"$BASE"-PBNC.txt\n- echo "${BASE},${BASE}-PBNC.txt" >> PBNC_file_paths.txt\n- \n- ##====== Sequence Duplication Level (SDL)\n- awk \'/^>>Sequence Duplication Levels/ {flag=1; next} /END_MODULE/ {flag=0} flag\' "$BASE"/fastqc_data.txt >"$BASE"-SDL.txt\n- echo "${BASE},${BASE}-SDL.txt" >> SDL_file_paths.txt\n- \n- ##====== Sequence Length Distribution (SLD)\n- awk \'/^>>Sequence Length Distribution/ {flag=1; next} /END_MODULE/ {flag=0} flag\' "$BASE"/fastqc_data.txt >"$BASE"-'..b')\n- PBNC_df = rbind(PBNC_df, pbnc_df)\n-}\n-```\n-\n-\n-```{r}\n-PBNC_df$N.Count = PBNC_df$N.Count * 100\n-max_phred = max(PBNC_df$N.Count) + 5\n-hchart(PBNC_df, "line", hcaes(x = as.character(Base), y = N.Count, group = sample_id)) %>%\n- hc_title(\n- text = "Per Base N Content"\n- ) %>%\n- hc_xAxis(\n- title = list(text = "Base Position")\n- ) %>%\n- hc_yAxis(\n- title = list(text = "N %"),\n- plotLines = list(\n- list(label = list(text = "N = 5%"),\n- width = 2,\n- dashStyle = "dash",\n- color = "red",\n- value = 5)\n- )\n- ) %>% \n- hc_exporting(enabled = TRUE)\n-```\n-\n-\n-\n-\n-## Per Sequence Quality Scores\n-\n-```{r}\n-PSQS_df = data.frame()\n-PSQS_file_paths = read.csv(\'PSQS_file_paths.txt\', \n- header = TRUE, stringsAsFactors = FALSE)\n-for(i in 1:nrow(PSQS_file_paths)) {\n- # file_path = paste0(\'REPORT_OUTPUT_DIR/\', PSQS_file_paths[i,2])\n- file_path = PSQS_file_paths[i,2]\n- psqs_df = read.csv(file_path,\n- sep=\'\\t\', header=TRUE, stringsAsFactors = FALSE) \n- psqs_df$sample_id = rep(PSQS_file_paths[i,1], nrow(psqs_df))\n- PSQS_df = rbind(PSQS_df, psqs_df)\n-}\n-```\n-\n-\n-```{r}\n-max_phred = max(PSQS_df$X.Quality) + 5\n-hchart(PSQS_df, "line", hcaes(x = X.Quality, y = Count, group = sample_id)) %>%\n- hc_title(\n- text = "Per Sequence Quality Score"\n- ) %>%\n- hc_xAxis(\n- title = list(text = "Mean Sequence Quality Score"),\n- min = 0,\n- max = max_phred,\n- plotLines = list(\n- list(label = list(text = "Phred Score = 27"),\n- width = 2,\n- dashStyle = "dash",\n- color = "green",\n- value = 27),\n- list(label = list(text = "Phred Score = 20"),\n- width = 2,\n- color = "red",\n- value = 20)\n- )\n- ) %>% \n- hc_exporting(enabled = TRUE)\n-```\n-\n-\n-## Per Sequence GC Content\n-\n-\n-```{r}\n-PSGC_df = data.frame()\n-PSGC_file_paths = read.csv(\'PSGC_file_paths.txt\', \n- header = TRUE, stringsAsFactors = FALSE)\n-for(i in 1:nrow(PSGC_file_paths)) {\n- # file_path = paste0(\'REPORT_OUTPUT_DIR/\', PSGC_file_paths[i,2])\n- file_path = PSGC_file_paths[i,2]\n- psgc_df = read.csv(file_path,\n- sep=\'\\t\', header=TRUE, stringsAsFactors = FALSE) \n- psgc_df$sample_id = rep(PSGC_file_paths[i,1], nrow(psgc_df))\n- PSGC_df = rbind(PSGC_df, psgc_df)\n-}\n-```\n-\n-\n-```{r}\n-max_phred = max(PSGC_df$Count) + 5\n-hchart(PSGC_df, "line", hcaes(x = X.GC.Content, y = Count, group = sample_id)) %>%\n- hc_title(\n- text = "Per Sequence GC Content"\n- ) %>%\n- hc_xAxis(\n- title = list(text = "% GC")\n- ) %>%\n- hc_exporting(enabled = TRUE)\n-```\n-\n-\n-## Per Base Sequence Content\n-\n-```{r}\n-PBSC_df = data.frame()\n-PBSC_file_paths = read.csv(\'PBSC_file_paths.txt\',\n- header = TRUE, stringsAsFactors = FALSE)\n-for(i in 1:nrow(PBSC_file_paths)) {\n- # file_path = paste0(\'REPORT_OUTPUT_DIR/\', PBSC_file_paths[i,2])\n- file_path = PBSC_file_paths[i,2]\n- pbsc_df = read.csv(file_path,\n- sep=\'\\t\', header=TRUE, stringsAsFactors = FALSE) %>%\n- mutate(Base1=as.numeric(str_split_fixed(X.Base, \'-\', 2)[,1]),\n- Base2=as.numeric(str_split_fixed(X.Base, \'-\', 2)[,2])) %>%\n- (function (df) {\n- df1 = select(df, -Base2)\n- df2 = select(df, -Base1) %>% filter(Base2 != \'\')\n- colnames(df1) = c(colnames(df1)[1:5], \'Base\')\n- colnames(df2) = c(colnames(df2)[1:5], \'Base\')\n- res = rbind(df1, df2) %>% arrange(Base)\n- return(res)\n- })\n- pbsc_df$sample_id = rep(PBSC_file_paths[i,1], nrow(pbsc_df))\n- PBSC_df = rbind(PBSC_df, pbsc_df)\n-}\n-```\n-\n-\n-```{r out.width="100%"}\n-PBSC_df_2 = select(PBSC_df, -X.Base) %>%\n- melt(id = c(\'Base\', \'sample_id\'), value.name = \'base_percentage\')\n-p = ggplot(data = PBSC_df_2, aes(x = Base, y = base_percentage, group = variable, color = variable)) +\n- geom_line() +\n- facet_wrap(~ sample_id)\n-ggplotly(p)\n-```\n-\n-\n-# Session Info\n-\n-```{r \'session info\'}\n-sessionInfo()\n-```\n-\n-\n' |
| b |
| diff -r 2efa46ce2c4c -r d1d20f341632 fastqc_report_render.R --- a/fastqc_report_render.R Wed Oct 18 22:06:39 2017 -0400 +++ b/fastqc_report_render.R Thu Oct 19 00:11:14 2017 -0400 |
| b |
| @@ -40,12 +40,13 @@ ##------- 1. input data --------------------- args_list$ECHO = c('echo', 'e', '1', 'character') args_list$READS = c('reads', 'r', '1', 'character') + args_list$NAMES = c('names', 'n', '1', 'character') ##--------2. output report and outputs -------------- - args_list$REPORT_HTML = c('report_html', 'r', '1', 'character') + args_list$REPORT_HTML = c('report_html', 'o', '1', 'character') args_list$REPORT_DIR = c('report_dir', 'd', '1', 'character') args_list$SINK_MESSAGE = c('sink_message', 's', '1', 'character') ##--------3. .Rmd templates in the tool directory ---------- - args_list$FASTQC_REPORT_RMD = c('fastqc_report_rmd', 't', '1', 'character') + args_list$FASTQC_REPORT_RMD = c('fastqc_report_rmd', 'p', '1', 'character') ##----------------------------------------------------------- opt = getopt(t(as.data.frame(args_list))) @@ -68,7 +69,7 @@ gsub('READS', opt$reads, x) }) %>% (function(x) { - gsub('REPORT_DIR', opt$output_dir, x) + gsub('REPORT_DIR', opt$report_dir, x) }) %>% (function(x) { fileConn = file('fastqc_report.Rmd') |
| b |
| diff -r 2efa46ce2c4c -r d1d20f341632 fastqc_report_render_ori.R --- a/fastqc_report_render_ori.R Wed Oct 18 22:06:39 2017 -0400 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 |
| [ |
| @@ -1,87 +0,0 @@ -##======= Handle arguments from command line ======== -# setup R error handline to go to stderr -options(show.error.messages = FALSE, -error = function(){ - cat(geterrmessage(), file = stderr()) - quit("no", 1, F) -}) - -# we need that to not crash galaxy with an UTF8 error on German LC settings. -loc = Sys.setlocale("LC_MESSAGES", "en_US.UTF-8") - -# suppress warning -options(warn = - 1) - -options(stringsAsFactors = FALSE, useFancyQuotes = FALSE) -args = commandArgs(trailingOnly = TRUE) - -suppressPackageStartupMessages({ - library(getopt) - library(tools) -}) - -# column 1: the long flag name -# column 2: the short flag alias. A SINGLE character string -# column 3: argument mask -# 0: no argument -# 1: argument required -# 2: argument is optional -# column 4: date type to which the flag's argument shall be cast. -# possible values: logical, integer, double, complex, character. -spec_list = list() -spec_list$READS = c('reads', 'r', '1', 'character') -spec_list$ECHO = c('echo', 'e', '1', 'character') -spec_list$FASTQC_TPL = c('fastqc_tpl', 'p', 1, 'character') -spec_list$REPORT = c('report', 'o', '1', 'character') -spec_list$REPORT_OUTPUT_DIR = c('report_output_dir', 'd', '1', 'character') - - -spec = t(as.data.frame(spec_list)) - -opt = getopt(spec) -# arguments are accessed by long flag name (the first column in the spec matrix) -# NOT by element name in the spec_list -# example: opt$help, opt$expression_file -##====== End of arguments handling ========== - - -mgsub = function(pattern, replacement, x) { - if (length(pattern) != length(replacement)) { - stop("pattern and replacement have to be the same in length") - } - - result = x - - for (i in 1 : length(pattern)) { - result = try(gsub(pattern[i], replacement[i], x = result)) - } - - result -} - - -##====== replace variables in tpl file ====== -p = c('READS', -'ECHO', -'FASTQC_TPL', -'REPORT_OUTPUT_DIR', -'REPORT') -r = c(opt$reads, -opt$echo, -opt$fastqc_tpl, -opt$report_output_dir, -opt$report) - -fastqc_report_tpl = mgsub(p, r, readLines(opt$fastqc_tpl)) - -##====== write replaced text into Rmd file === -fileConn = file('fastqc_report.Rmd') -writeLines(fastqc_report_tpl, con = fileConn) -close(fileConn) - -##====== render Rmd files ==================== -rmarkdown::render('fastqc_report.Rmd') -file.copy('fastqc_report.html', opt$report, recursive = TRUE) -paste0('cp -r ./* ', opt$report_output_dir) %>% -system() - |