Mercurial > repos > mingchen0919 > aurora_deseq2_site

changeset 0:6f94b4b9de44 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload
author mingchen0919
date Tue, 27 Feb 2018 23:57:53 -0500
parents
children 32210899a3dd
files DESeq.Rmd DESeq.xml DESeq_01.Rmd DESeq_02.Rmd DESeq_03.Rmd DESeq_04.Rmd DESeq_05.Rmd DESeq_index.Rmd DESeq_render.R DESeq_results.Rmd DESeq_results.xml DESeq_results_01.Rmd DESeq_results_02.Rmd DESeq_results_03.Rmd DESeq_results_04.Rmd DESeq_results_index.Rmd DESeq_results_render.R DESeq_results_site.yml DESeq_site.yml
diffstat 17 files changed, 887 insertions(+), 0 deletions(-) [+]
line wrap: on
line diff
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/DESeq.Rmd	Tue Feb 27 23:57:53 2018 -0500
@@ -0,0 +1,97 @@
+---
+title: 'DESeq2: Perform DESeq analysis'
+output:
+    html_document:
+      number_sections: true
+      toc: true
+      theme: cosmo
+      highlight: tango
+---
+
+```{r setup, include=FALSE, warning=FALSE, message=FALSE}
+knitr::opts_chunk$set(
+  echo = as.logical(opt$X_e),
+  error = TRUE
+)
+```
+
+# `DESeqDataSet` object
+
+```{r 'DESeqDataSet object'}
+count_file_paths = strsplit(opt$X_P, ',')[[1]]
+count_file_names = strsplit(opt$X_N, ',')[[1]]
+sample_table = read.table(opt$X_S, header = TRUE)
+row.names(sample_table) = sample_table[,2]
+sample_table = sample_table[count_file_names, ]
+
+## copy count files into OUTPUT_DIR/counts
+dir.create(paste0(OUTPUT_DIR, '/counts'), recursive = TRUE)
+file_copy = file.copy(count_file_paths, paste0(OUTPUT_DIR, '/counts/', count_file_names), overwrite = TRUE)
+
+## DESeqDataSet object
+dds = DESeqDataSetFromHTSeqCount(sampleTable = sample_table,
+                                 directory = paste0(OUTPUT_DIR, '/counts'),
+                                 design = formula(opt$X_p))
+dds
+```
+
+# Pre-filtering the dataset.
+
+We can remove the rows that have 0 or 1 count to reduce object size and increase the calculation speed.
+
+* Number of rows before pre-filtering
+```{r}
+nrow(dds)
+```
+
+* Number of rows after pre-filtering
+```{r}
+dds = dds[rowSums(counts(dds)) > 1, ]
+nrow(dds)
+```
+
+# Peek at data {.tabset}
+
+## Count Data
+
+```{r 'count data'}
+datatable(head(counts(dds), 100), style="bootstrap", 
+          class="table-condensed", options = list(dom = 'tp', scrollX = TRUE))
+```
+
+## Sample Table 
+
+```{r 'sample table'}
+datatable(sample_table, style="bootstrap",
+          class="table-condensed", options = list(dom = 'tp', scrollX = TRUE))
+```
+
+# Sample distance on variance stabilized data {.tabset}
+
+## `rlog` Stabilizing transformation
+
+```{r}
+rld = rlog(dds, blind = FALSE)
+datatable(head(assay(rld), 100), style="bootstrap", 
+          class="table-condensed", options = list(dom = 'tp', scrollX = TRUE))
+```
+
+## Sample distance
+
+```{r}
+sampleDists <- dist(t(assay(rld)))
+sampleDists
+```
+
+# Differential expression analysis
+
+```{r}
+dds <- DESeq(dds)
+```
+
+```{r echo=FALSE}
+# save objects except for opt.
+save(list=ls()[ls() != "opt"], file=opt$X_w)
+```
+
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/DESeq.xml	Tue Feb 27 23:57:53 2018 -0500
@@ -0,0 +1,106 @@
+<tool name="DESeq2: Analysis" id='deseq2' version="2.0.1">
+    <description>
+        perform differential expression analysis
+    </description>
+    <requirements>
+        <requirement type="package" version="1.15.0.6-0">pandoc</requirement>
+        <requirement type="package" version="1.20.0">r-getopt</requirement>
+        <requirement type="package" version="1.6">r-rmarkdown</requirement>
+        <requirement type="package" version="1.14.1">bioconductor-deseq2</requirement>
+        <requirement type="package" version="0.2">r-dt</requirement>
+        <requirement type="package" version="1.0.8">r-pheatmap</requirement>
+    </requirements>
+    <stdio>
+        <regex match="XXX" source="stderr" level="warning"
+               description="Check the warnings_and_errors.txt file for more details."/>
+    </stdio>
+    <command><![CDATA[
+        
+
+        Rscript '${__tool_directory__}/DESeq_render.R'
+
+            -e $echo
+
+		    -o $report
+		    -d $report.files_path
+		    -s $sink_message
+		    -t '${__tool_directory__}'
+
+
+		    ##----- code chunk to get file paths and raw file names for a multiple inputs data field ----
+            #set $sep = ''
+            #set $count_file_paths = ''
+            #set $count_file_names = ''
+            #for $count_file in $count_files:
+                #set $count_file_paths += $sep + str($count_file)
+                #set $count_file_names += $sep + str($count_file.name)
+                #set $sep = ','
+            #end for
+            ##----------------- end for getting file names and file paths ------------------------------
+            -P '$count_file_paths'
+            -N '$count_file_names'
+            -S $sample_table
+            -p '$design_formula'
+            -w $deseq_workspace
+
+        
+    ]]></command>
+    <inputs>
+        <param type="boolean" name="echo" truevalue="TRUE" falsevalue="FALSE" checked="false"
+               label="Display analysis code in report?"/>
+        <param type="data" name="count_files" format="txt" multiple="true" label="Count files from htseq-count"/>
+        <param type="data" name="sample_table" format="txt" multiple="false" label="sample table file"
+               help="The sample table file contains a table. The first column is the sample name, the second column is
+                    the count file name and the rest of columns are treatment columns. The file names in this table have
+                    to be in the same order as the count files uploaded in the previous step. "/>
+        <param type="text" name="design_formula" value="~ condition_1 + condition_2" label="Design formula"
+               help="The simplest design formula for differential expression would be ~ condition, where condition
+                     is a column in colData(dds) that specifies which of two (or more groups) the samples belong to">
+            <sanitizer>
+                <valid initial="default">
+                    <add preset="string.printable"/>
+                    <add value="~"/>
+                </valid>
+            </sanitizer>
+        </param>
+    </inputs>
+    <outputs>
+        <data name="report" format="html" label="DESeq Analysis on ${on_string}"/>
+        <data format="txt" name="sink_message" label="Warnings and Errors on"
+              from_work_dir="warnings_and_errors.txt"/>
+        <data name="deseq_workspace" format="rdata" label="R workspace: DESeq analysis on ${on_string}"/>
+    </outputs>
+    <citations>
+        <citation type="bibtex">
+            @article{love2014moderated,
+            title={Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2},
+            author={Love, Michael I and Huber, Wolfgang and Anders, Simon},
+            journal={Genome biology},
+            volume={15},
+            number={12},
+            pages={550},
+            year={2014},
+            publisher={BioMed Central}
+            }
+        </citation>
+        <citation type="bibtex"><![CDATA[
+            @article{allaire2016rmarkdown,
+            title={rmarkdown: Dynamic Documents for R, 2016},
+            author={Allaire, J and Cheng, Joe and Xie, Yihui and McPherson, Jonathan and Chang, Winston and Allen, Jeff
+            and Wickham, Hadley and Atkins, Aron and Hyndman, Rob},
+            journal={R package version 0.9},
+            volume={6},
+            year={2016}
+            }
+        ]]></citation>
+        <citation type="bibtex"><![CDATA[
+            @book{xie2015dynamic,
+            title={Dynamic Documents with R and knitr},
+            author={Xie, Yihui},
+            volume={29},
+            year={2015},
+            publisher={CRC Press}
+            }
+        ]]></citation>
+    </citations>
+</tool>
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/DESeq_01.Rmd	Tue Feb 27 23:57:53 2018 -0500
@@ -0,0 +1,30 @@
+---
+output: html_document
+---
+
+```{r setup, include=FALSE, warning=FALSE, message=FALSE}
+knitr::opts_chunk$set(
+  echo = as.logical(opt$X_e),
+  error = TRUE
+)
+```
+
+# `DESeqDataSet` object
+
+```{r 'DESeqDataSet object'}
+count_file_paths = strsplit(opt$X_P, ',')[[1]]
+count_file_names = strsplit(opt$X_N, ',')[[1]]
+sample_table = read.table(opt$X_S, header = TRUE)
+row.names(sample_table) = sample_table[,2]
+sample_table = sample_table[count_file_names, ]
+
+## copy count files into OUTPUT_DIR/counts
+dir.create(paste0(OUTPUT_DIR, '/counts'), recursive = TRUE)
+file_copy = file.copy(count_file_paths, paste0(OUTPUT_DIR, '/counts/', count_file_names), overwrite = TRUE)
+
+## DESeqDataSet object
+dds = DESeqDataSetFromHTSeqCount(sampleTable = sample_table,
+                                 directory = paste0(OUTPUT_DIR, '/counts'),
+                                 design = formula(opt$X_p))
+dds
+```
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/DESeq_02.Rmd	Tue Feb 27 23:57:53 2018 -0500
@@ -0,0 +1,27 @@
+---
+output: html_document
+---
+
+```{r setup, include=FALSE, warning=FALSE, message=FALSE}
+knitr::opts_chunk$set(
+  echo = as.logical(opt$X_e),
+  error = TRUE
+)
+```
+
+
+
+# Pre-filtering the dataset.
+
+We can remove the rows that have 0 or 1 count to reduce object size and increase the calculation speed.
+
+* Number of rows before pre-filtering
+```{r}
+nrow(dds)
+```
+
+* Number of rows after pre-filtering
+```{r}
+dds = dds[rowSums(counts(dds)) > 1, ]
+nrow(dds)
+```
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/DESeq_03.Rmd	Tue Feb 27 23:57:53 2018 -0500
@@ -0,0 +1,27 @@
+---
+output: html_document
+---
+
+```{r setup, include=FALSE, warning=FALSE, message=FALSE}
+knitr::opts_chunk$set(
+  echo = as.logical(opt$X_e),
+  error = TRUE
+)
+```
+
+
+# Peek at data {.tabset}
+
+## Count Data
+
+```{r 'count data'}
+datatable(head(counts(dds), 100), style="bootstrap", 
+          class="table-condensed", options = list(dom = 'tp', scrollX = TRUE))
+```
+
+## Sample Table 
+
+```{r 'sample table'}
+datatable(sample_table, style="bootstrap",
+          class="table-condensed", options = list(dom = 'tp', scrollX = TRUE))
+```
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/DESeq_04.Rmd	Tue Feb 27 23:57:53 2018 -0500
@@ -0,0 +1,29 @@
+---
+output: html_document
+---
+
+```{r setup, include=FALSE, warning=FALSE, message=FALSE}
+knitr::opts_chunk$set(
+  echo = as.logical(opt$X_e),
+  error = TRUE
+)
+```
+
+
+
+# Sample distance on variance stabilized data {.tabset}
+
+## `rlog` Stabilizing transformation
+
+```{r}
+rld = rlog(dds, blind = FALSE)
+datatable(head(assay(rld), 100), style="bootstrap", 
+          class="table-condensed", options = list(dom = 'tp', scrollX = TRUE))
+```
+
+## Sample distance
+
+```{r}
+sampleDists <- dist(t(assay(rld)))
+sampleDists
+```
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/DESeq_05.Rmd	Tue Feb 27 23:57:53 2018 -0500
@@ -0,0 +1,23 @@
+---
+output: html_document
+---
+
+```{r setup, include=FALSE, warning=FALSE, message=FALSE}
+knitr::opts_chunk$set(
+  echo = as.logical(opt$X_e),
+  error = TRUE
+)
+```
+
+
+# Differential expression analysis
+
+```{r}
+dds <- DESeq(dds)
+```
+
+```{r echo=FALSE}
+# save useful objects.
+save(dds, rld, sample_table, sampleDists, zz, file=opt$X_w)
+```
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/DESeq_render.R	Tue Feb 27 23:57:53 2018 -0500
@@ -0,0 +1,86 @@
+##============ Sink warnings and errors to a file ==============
+## use the sink() function to wrap all code within it.
+##==============================================================
+zz = file('warnings_and_errors.txt')
+sink(zz)
+sink(zz, type = 'message')
+
+#------------import libraries--------------------
+options(stringsAsFactors = FALSE)
+
+library(getopt)
+library(rmarkdown)
+library(DESeq2)
+library(pheatmap)
+library(DT)
+library(ggplot2)
+library(genefilter)
+library(RColorBrewer)
+#------------------------------------------------
+
+
+#------------get arguments into R--------------------
+# getopt_specification_matrix(extract_short_flags('fastqc_report.xml')) %>%
+#   write.table(file = 'spec.txt', sep = ',', row.names = FALSE, col.names = TRUE, quote = FALSE)
+
+
+spec_matrix = as.matrix(
+  data.frame(stringsAsFactors=FALSE,
+             long_flags = c("X_e", "X_o", "X_d", "X_s", "X_t", "X_P", "X_N",
+                            "X_S", "X_p", "X_w"),
+             short_flags = c("e", "o", "d", "s", "t", "P", "N", "S", "p", "w"),
+             argument_mask_flags = c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L),
+             data_type_flags = c("character", "character", "character", "character",
+                                 "character", "character", "character",
+                                 "character", "character", "character")
+  )
+)
+opt = getopt(spec_matrix)
+#----------------------------------------------------
+
+
+#-----------using passed arguments in R 
+#           to define system environment variables---
+do.call(Sys.setenv, opt[-1])
+#----------------------------------------------------
+
+#---------- often used variables ----------------
+# OUTPUT_REPORT: path to galaxy output report
+# OUTPUT_DIR: path to the output associated directory, which stores all outputs
+# TOOL_DIR: path to the tool installation directory
+OUTPUT_DIR = opt$X_d
+TOOL_DIR =   opt$X_t
+OUTPUT_REPORT = opt$X_o
+
+
+# create the output associated directory to store all outputs
+dir.create(OUTPUT_DIR, recursive = TRUE)
+
+#-----------------render site--------------
+# copy site generating materials into OUTPUT_DIR
+dir.create(paste0(OUTPUT_DIR, '/site_generator'), recursive = TRUE)
+command_cp = paste0('cp -r ', TOOL_DIR, '/DESeq_0*.Rmd ', OUTPUT_DIR, '/site_generator')
+system(command_cp)
+system(paste0('cp -r ', TOOL_DIR, '/DESeq_site.yml ', OUTPUT_DIR, '/site_generator/_site.yml'))
+system(paste0('cp -r ', TOOL_DIR, '/DESeq_index.Rmd ', OUTPUT_DIR, '/site_generator/index.Rmd'))
+# render site to OUTPUT_DIR/_site, this is configured in the "_site.yml" file
+dir.create(paste0(OUTPUT_DIR, '/_site'))
+render_site(input = paste0(OUTPUT_DIR, '/site_generator'))
+# remove site generating materials from output associated directory
+print(unlink(paste0(OUTPUT_DIR, '/site_generator'), recursive = TRUE))
+# move _site/* into output associated directory
+move_cmd = paste0('mv ', OUTPUT_DIR, '/_site/* ', OUTPUT_DIR)
+system(move_cmd)
+#------------------------------------------
+
+#-----link index.html to output-----
+cp_index = paste0('cp ', OUTPUT_DIR, '/index.html ', OUTPUT_REPORT)
+system(cp_index)
+#-----------------------------------
+
+#==============the end==============
+
+
+##--------end of code rendering .Rmd templates----------------
+sink()
+##=========== End of sinking output=============================
\ No newline at end of file
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/DESeq_results.Rmd	Tue Feb 27 23:57:53 2018 -0500
@@ -0,0 +1,109 @@
+---
+title: 'DESeq2: Results'
+output:
+    html_document:
+      number_sections: true
+      toc: true
+      theme: cosmo
+      highlight: tango
+---
+
+```{r setup, include=FALSE, warning=FALSE, message=FALSE}
+knitr::opts_chunk$set(
+  echo = as.logical(opt$X_e),
+  error = TRUE
+)
+```
+
+
+```{r eval=TRUE}
+# Import workspace
+# fcp = file.copy(opt$X_W, "deseq.RData")
+load(opt$X_W)
+```
+
+# Results {.tabset}
+
+## Result table
+
+```{r}
+cat('--- View the top 100 rows of the result table ---')
+res <- results(dds, contrast = c(opt$X_C, opt$X_T, opt$X_K))
+write.csv(as.data.frame(res), file = opt$X_R)
+res_df = as.data.frame(res)[1:100, ]
+datatable(res_df, style="bootstrap", filter = 'top',
+          class="table-condensed", options = list(dom = 'tp', scrollX = TRUE))
+```
+
+## Result summary
+
+```{r}
+summary(res)
+```
+
+
+# MA-plot {.tabset}
+
+
+
+```{r}
+cat('--- Shrinked with Bayesian procedure ---')
+plotMA(res)
+```
+
+
+# Histogram of p values
+
+```{r}
+hist(res$pvalue[res$baseMean > 1], breaks = 0:20/20,
+     col = "grey50", border = "white", main = "",
+     xlab = "Mean normalized count larger than 1")
+```
+
+
+# Visualization {.tabset}
+## Gene clustering
+
+```{r}
+clustering_groups = strsplit(opt$X_M, ',')[[1]]
+
+topVarGenes <- head(order(rowVars(assay(rld)), decreasing = TRUE), 20)
+mat  <- assay(rld)[ topVarGenes, ]
+mat  <- mat - rowMeans(mat)
+annotation_col <- as.data.frame(colData(rld)[, clustering_groups])
+colnames(annotation_col) = clustering_groups
+rownames(annotation_col) = colnames(mat)
+pheatmap(mat, annotation_col = annotation_col)
+```
+
+## Sample-to-sample distance
+
+```{r}
+sampleDistMatrix <- as.matrix( sampleDists )
+colors <- colorRampPalette( rev(brewer.pal(9, "Blues")) )(255)
+pheatmap(sampleDistMatrix,
+         clustering_distance_cols = sampleDists,
+         col = colors)
+```
+
+## PCA plot 
+
+```{r}
+plotPCA(rld, intgroup = clustering_groups)
+```
+
+## MDS plot {.tabset}
+
+### Data table
+```{r}
+mds <- as.data.frame(colData(rld))  %>%
+         cbind(cmdscale(sampleDistMatrix))
+knitr::kable(mds)
+```
+
+### Plot
+```{r}
+ggplot(mds, aes(x = `1`, y = `2`, col = time)) +
+  geom_point(size = 3) + coord_fixed()
+```
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/DESeq_results.xml	Tue Feb 27 23:57:53 2018 -0500
@@ -0,0 +1,99 @@
+<tool id="deseq_results" name="DESeq2: Results" version="2.0.1">
+    <requirements>
+        <requirement type="package" version="1.15.0.6-0">pandoc</requirement>
+        <requirement type="package" version="1.20.0">r-getopt</requirement>
+        <requirement type="package" version="1.6">r-rmarkdown</requirement>
+        <requirement type="package" version="1.14.1">bioconductor-deseq2</requirement>
+        <requirement type="package" version="0.2">r-dt</requirement>
+        <requirement type="package" version="1.0.8">r-pheatmap</requirement>
+    </requirements>
+    <description>
+        display DESeq2 analysis results.
+    </description>
+    <stdio>
+        <!--redirecting stderr to a file. "XXX" is used to match with nothing so that tool running won't be interrupted during testing-->
+        <regex match="XXX"
+               source="stderr"
+               level="warning"
+               description="Check the warnings_and_errors.txt file for more details."/>
+    </stdio>
+    <command>
+        <![CDATA[
+
+        Rscript '${__tool_directory__}/DESeq_results_render.R'
+
+            -e $echo
+		    -o $report
+		    -d $report.files_path
+		    -s $sink_message
+		    -t '${__tool_directory__}'
+
+            ## 1. input data
+
+            -W $deseq_workspace
+            -C '$contrast_factor'
+            -T '$treatment'
+            -K '$condition'
+
+            -M '$clustering_factors'
+
+            ## 2. output report and report site directory
+            -R $deseq_results
+
+
+        ]]>
+    </command>
+    <inputs>
+        <param type="boolean" name="echo" truevalue="TRUE" falsevalue="FALSE" checked="false"
+               label="Display analysis code in report?"/>
+        <param type="data" name="deseq_workspace" format="rdata" multiple="false" optional="false"
+               label="Workspace from tool DESeq2: DESeq"/>
+        <param type="text" name="contrast_factor" label="Factor" optional="false"
+               help="the name of a factor in the design formula"/>
+        <param type="text" name="treatment" label="Treatment level" optional="false"
+               help=" the name of the numerator level for the fold change"/>
+        <param type="text" name="condition" label="Condition level" optional="false"
+               help=" the name of the denominator level for the fold change"/>
+        <param type="text" name="clustering_factors" title="Gene clustering factors" optional="false"
+               label="factors of interest for clustering samples and PCA plot"
+               help="A single factor or multiple factors from the design formula. Multiple factors are separated by comma (,)."/>
+    </inputs>
+    <outputs>
+        <data format="html" name="report" label="DESeq results report on ${on_string}" />
+        <data format="txt" name="sink_message" label="Warnings and Errors" from_work_dir="warnings_and_errors.txt"/>
+        <data format="csv" name="deseq_results" label="DESeq results on ${on_string}" from_work_dir="deseq_results.csv" />
+    </outputs>
+    <citations>
+        <citation type="bibtex">
+            @article{love2014moderated,
+            title={Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2},
+            author={Love, Michael I and Huber, Wolfgang and Anders, Simon},
+            journal={Genome biology},
+            volume={15},
+            number={12},
+            pages={550},
+            year={2014},
+            publisher={BioMed Central}
+            }
+        </citation>
+        <citation type="bibtex">
+            @article{allaire2016rmarkdown,
+            title={rmarkdown: Dynamic Documents for R, 2016},
+            author={Allaire, J and Cheng, Joe and Xie, Yihui and McPherson, Jonathan and Chang, Winston and Allen, Jeff
+            and Wickham, Hadley and Atkins, Aron and Hyndman, Rob},
+            journal={R package version 0.9},
+            volume={6},
+            year={2016}
+            }
+        </citation>
+        <citation type="bibtex">
+            @book{xie2015dynamic,
+            title={Dynamic Documents with R and knitr},
+            author={Xie, Yihui},
+            volume={29},
+            year={2015},
+            publisher={CRC Press}
+            }
+        </citation>
+    </citations>
+</tool>
\ No newline at end of file
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/DESeq_results_01.Rmd	Tue Feb 27 23:57:53 2018 -0500
@@ -0,0 +1,30 @@
+---
+output: html_document
+---
+
+```{r setup, include=FALSE, warning=FALSE, message=FALSE}
+knitr::opts_chunk$set(
+  echo = as.logical(opt$X_e),
+  error = TRUE
+)
+```
+
+
+```{r eval=TRUE}
+# Import workspace
+# fcp = file.copy(opt$X_W, "deseq.RData")
+load(opt$X_W)
+```
+
+# Results {.tabset}
+
+## Result table
+
+```{r}
+cat('--- View the top 100 rows of the result table ---')
+res <- results(dds, contrast = c(opt$X_C, opt$X_T, opt$X_K))
+write.csv(as.data.frame(res), file = opt$X_R)
+res_df = as.data.frame(res)[1:100, ]
+datatable(res_df, style="bootstrap", filter = 'top',
+          class="table-condensed", options = list(dom = 'tp', scrollX = TRUE))
+```
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/DESeq_results_02.Rmd	Tue Feb 27 23:57:53 2018 -0500
@@ -0,0 +1,17 @@
+---
+output: html_document
+---
+
+```{r setup, include=FALSE, warning=FALSE, message=FALSE}
+knitr::opts_chunk$set(
+  echo = as.logical(opt$X_e),
+  error = TRUE
+)
+```
+
+# MA-plot {.tabset}
+
+```{r}
+cat('--- Shrinked with Bayesian procedure ---')
+plotMA(res)
+```
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/DESeq_results_03.Rmd	Tue Feb 27 23:57:53 2018 -0500
@@ -0,0 +1,20 @@
+---
+output: html_document
+---
+
+```{r setup, include=FALSE, warning=FALSE, message=FALSE}
+knitr::opts_chunk$set(
+  echo = as.logical(opt$X_e),
+  error = TRUE
+)
+```
+
+
+# Histogram of p values
+
+```{r}
+hist(res$pvalue[res$baseMean > 1], breaks = 0:20/20,
+     col = "grey50", border = "white", main = "",
+     xlab = "Mean normalized count larger than 1")
+```
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/DESeq_results_04.Rmd	Tue Feb 27 23:57:53 2018 -0500
@@ -0,0 +1,57 @@
+---
+output: html_document
+---
+
+```{r setup, include=FALSE, warning=FALSE, message=FALSE}
+knitr::opts_chunk$set(
+  echo = as.logical(opt$X_e),
+  error = TRUE
+)
+```
+
+
+# Visualization {.tabset}
+## Gene clustering
+
+```{r}
+clustering_groups = strsplit(opt$X_M, ',')[[1]]
+
+topVarGenes <- head(order(rowVars(assay(rld)), decreasing = TRUE), 20)
+mat  <- assay(rld)[ topVarGenes, ]
+mat  <- mat - rowMeans(mat)
+annotation_col <- as.data.frame(colData(rld)[, clustering_groups])
+colnames(annotation_col) = clustering_groups
+rownames(annotation_col) = colnames(mat)
+pheatmap(mat, annotation_col = annotation_col)
+```
+
+## Sample-to-sample distance
+
+```{r}
+sampleDistMatrix <- as.matrix( sampleDists )
+colors <- colorRampPalette( rev(brewer.pal(9, "Blues")) )(255)
+pheatmap(sampleDistMatrix,
+         clustering_distance_cols = sampleDists,
+         col = colors)
+```
+
+## PCA plot 
+
+```{r}
+plotPCA(rld, intgroup = clustering_groups)
+```
+
+## MDS plot {.tabset}
+
+### Data table
+```{r}
+mds <- as.data.frame(colData(rld))  %>%
+         cbind(cmdscale(sampleDistMatrix))
+knitr::kable(mds)
+```
+
+### Plot
+```{r}
+ggplot(mds, aes(x = `1`, y = `2`, col = time)) +
+  geom_point(size = 3) + coord_fixed()
+```
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/DESeq_results_render.R	Tue Feb 27 23:57:53 2018 -0500
@@ -0,0 +1,88 @@
+##============ Sink warnings and errors to a file ==============
+## use the sink() function to wrap all code within it.
+##==============================================================
+zz = file('warnings_and_errors.txt')
+sink(zz)
+sink(zz, type = 'message')
+
+#------------import libraries--------------------
+options(stringsAsFactors = FALSE)
+
+library(getopt)
+library(rmarkdown)
+library(DESeq2)
+library(pheatmap)
+library(DT)
+library(ggplot2)
+library(genefilter)
+library(RColorBrewer)
+#------------------------------------------------
+
+
+#------------get arguments into R--------------------
+# getopt_specification_matrix(extract_short_flags('fastqc_report.xml')) %>%
+#   write.table(file = 'spec.txt', sep = ',', row.names = FALSE, col.names = TRUE, quote = FALSE)
+
+
+spec_matrix = as.matrix(
+  data.frame(stringsAsFactors=FALSE,
+             long_flags = c("X_e", "X_W", "X_C", "X_T", "X_K", "X_M", "X_o",
+                            "X_d", "X_s", "X_R", "X_t"),
+             short_flags = c("e", "W", "C", "T", "K", "M", "o", "d", "s", "R",
+                             "t"),
+             argument_mask_flags = c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L),
+             data_type_flags = c("character", "character", "character", "character",
+                                 "character", "character", "character",
+                                 "character", "character", "character", "character")
+  )
+)
+opt = getopt(spec_matrix)
+opt
+#----------------------------------------------------
+
+
+#-----------using passed arguments in R 
+#           to define system environment variables---
+do.call(Sys.setenv, opt[-1])
+#----------------------------------------------------
+
+#---------- often used variables ----------------
+# OUTPUT_REPORT: path to galaxy output report
+# OUTPUT_DIR: path to the output associated directory, which stores all outputs
+# TOOL_DIR: path to the tool installation directory
+OUTPUT_DIR = opt$X_d
+TOOL_DIR =   opt$X_t
+OUTPUT_REPORT = opt$X_o
+
+
+# create the output associated directory to store all outputs
+dir.create(OUTPUT_DIR, recursive = TRUE)
+
+#-----------------render site--------------
+# copy site generating materials into OUTPUT_DIR
+dir.create(paste0(OUTPUT_DIR, '/site_generator'), recursive = TRUE)
+command_cp = paste0('cp -r ', TOOL_DIR, '/DESeq_results_*.Rmd ', OUTPUT_DIR, '/site_generator')
+system(command_cp)
+system(paste0('cp -r ', TOOL_DIR, '/DESeq_results_site.yml ', OUTPUT_DIR, '/site_generator/_site.yml'))
+system(paste0('cp -r ', TOOL_DIR, '/DESeq_results_index.Rmd ', OUTPUT_DIR, '/site_generator/index.Rmd'))
+# render site to OUTPUT_DIR/_site, this is configured in the "_site.yml" file
+dir.create(paste0(OUTPUT_DIR, '/_site'))
+render_site(input = paste0(OUTPUT_DIR, '/site_generator'))
+# remove site generating materials from output associated directory
+print(unlink(paste0(OUTPUT_DIR, '/site_generator'), recursive = TRUE))
+# move _site/* into output associated directory
+move_cmd = paste0('mv ', OUTPUT_DIR, '/_site/* ', OUTPUT_DIR)
+system(move_cmd)
+#------------------------------------------
+
+#-----link index.html to output-----
+cp_index = paste0('cp ', OUTPUT_DIR, '/index.html ', OUTPUT_REPORT)
+system(cp_index)
+#-----------------------------------
+
+#==============the end==============
+
+
+##--------end of code rendering .Rmd templates----------------
+sink()
+##=========== End of sinking output=============================
\ No newline at end of file
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/DESeq_results_site.yml	Tue Feb 27 23:57:53 2018 -0500
@@ -0,0 +1,21 @@
+name: "Analysis Report"
+output_dir: "../_site"
+navbar:
+    title: ""
+    type: inverse
+    left:
+        - text: "Home"
+          icon: fa-home
+          href: index.html
+        - text: "Results"
+          href: DESeq_results_01.html
+        - text: "MA-plot"
+          href: DESeq_results_02.html
+        - text: "Histogram of p values"
+          href: DESeq_results_03.html
+        - text: "Visualization"
+          href: DESeq_results_04.html
+output:
+  html_document:
+    theme: cosmo
+    highlight: textmate
\ No newline at end of file
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/DESeq_site.yml	Tue Feb 27 23:57:53 2018 -0500
@@ -0,0 +1,21 @@
+name: "Analysis Report"
+output_dir: "../_site"
+navbar:
+    title: ""
+    type: inverse
+    left:
+        - text: "Home"
+          icon: fa-home
+          href: index.html
+        - text: "DESeqDataSet"
+          href: DESeq_01.html
+        - text: "Pre-filtering"
+          href: DESeq_02.html
+        - text: "Data"
+          href: DESeq_03.html
+        - text: "Sample distance"
+          href: DESeq_04.html
+output:
+  html_document:
+    theme: cosmo
+    highlight: textmate
\ No newline at end of file