--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/rmarkdown_report.Rmd	Fri Dec 14 00:21:26 2018 -0500
@@ -0,0 +1,156 @@
+---
+title: 'Aurora DESeq2 Report'
+output:
+    html_document:
+      highlight: pygments
+---
+
+```{r setup, include=FALSE, warning=FALSE, message=FALSE}
+knitr::opts_chunk$set(error = TRUE, echo = FALSE)
+```
+
+```{css, echo=FALSE}
+pre code, pre, code {
+  white-space: pre !important;
+  overflow-x: scroll !important;
+  word-break: keep-all !important;
+  word-wrap: initial !important;
+}
+```
+
+```{r, echo=FALSE}
+# to make the css theme to work, <link></link> tags cannot be added directly 
+# as <script></script> tags as below.
+# it has to be added using a code chunk with the htmltool functions!!!
+css_link = tags$link()
+css_link$attribs = list(rel="stylesheet", href="vakata-jstree-3.3.5/dist/themes/default/style.min.css")
+css_link
+```
+
+```{r, eval=FALSE, echo=FALSE}
+# this code chunk is purely for adding comments
+# below is to add jQuery and jstree javascripts
+```
+
+<script src="vakata-jstree-3.3.5/dist/jstree.min.js"></script>
+
+
+```{r, eval=FALSE, echo=FALSE}
+# this code chunk is purely for adding comments
+# javascript code below is to build the file tree interface
+# see this for how to implement opening hyperlink: https://stackoverflow.com/questions/18611317/how-to-get-i-get-leaf-nodes-in-jstree-to-open-their-hyperlink-when-clicked-when
+```
+<script>
+  jQuery(function () {
+    // create an instance when the DOM is ready
+    jQuery('#jstree').jstree().bind("select_node.jstree", function (e, data) {
+     window.open( data.node.a_attr.href, data.node.a_attr.target )
+    });
+  });
+</script>
+
+
+```{r, eval=FALSE, echo=FALSE}
+---
+# ADD YOUR DATA ANALYSIS CODE AND MARKUP TEXT BELOW TO EXTEND THIS R MARKDOWN FILE
+---
+```
+
+## DESeq2 analysis
+
+```{r echo=FALSE}
+# import count data
+count_data = read.csv(opt$X_A, row.names = 1, header = TRUE)
+# import column data
+coldata = read.csv(opt$X_B, row.names = 1, header = TRUE)[colnames(count_data),,drop=FALSE]
+```
+
+```{r}
+f = gsub('~', '~ 1 +', opt$X_C) # build formula
+dds = DESeqDataSetFromMatrix(countData = count_data,
+                             colData = coldata,
+                             design = formula(f))
+
+
+# prefiltering
+keep <- rowSums(counts(dds)) >= 10
+dds <- dds[keep,]
+
+# Run DESeq
+if (opt$X_T == 'LRT') {
+  reduced_f = gsub(paste0('\\+\\s*', opt$X_D), '', f)
+  dds =  DESeq(dds, test=opt$X_T, fitType = opt$X_H, reduced = formula(reduced_f))
+} else {
+  dds = DESeq(dds, test=opt$X_T, fitType = opt$X_H)
+}
+
+## Differential expression test results
+res = results(dds, contrast = c(opt$X_D, opt$X_E, opt$X_F), alpha = opt$X_I)
+DT::datatable(as.data.frame(res))
+```
+
+
+```{r}
+# save all padj sorted res to tool output directory
+padj_sorted_res = res[order(res$padj), ]
+write.table(padj_sorted_res,
+            file = paste0(opt$X_d, '/padj-sorted-genes.txt'),
+            quote = FALSE)
+
+# save significant genes to a file in tool output directory
+sig_res = res[(res$padj < opt$X_I) & !is.na(res$padj), ]
+sig_res_sorted = sig_res[order(sig_res$padj), ]
+sig_res_sorted$feature_id = rownames(sig_res_sorted)
+n_col = ncol(sig_res_sorted)
+sig_res_sorted = sig_res_sorted[, c(n_col, 1:(n_col - 1))]
+write.table(sig_res_sorted, 
+            file = paste0(opt$X_d, '/padj-sorted-significant-genes.txt'), 
+            quote = FALSE, row.names = FALSE)
+```
+
+## MA-plot
+
+```{r warning=FALSE}
+log_fold_change = res$log2FoldChange
+base_mean = res$baseMean
+significant = res$padj
+significant[significant < 0.1] = 'yes'
+significant[significant != 'yes'] = 'no'
+
+maplot_df = data.frame(log_fold_change, base_mean, significant)
+maplot_df = maplot_df[!is.na(maplot_df$significant), ]
+p = ggplot(data = maplot_df) +
+  geom_point(mapping = aes(log(base_mean), log_fold_change, color = significant),
+             size = 0.5) +
+  scale_color_manual(name = 'Significant',
+                     values = c('no' = 'black', 'yes' = 'red'),
+                     labels = c('No', 'Yes')) +
+  xlab('Log base mean') +
+  ylab('Log fold change') +
+  theme_classic()
+
+plotly::ggplotly(p)
+```
+
+
+## Heatmap of count matrix
+
+```{r}
+ntd <- normTransform(dds)
+select <- order(rowMeans(counts(dds,normalized=TRUE)),
+                decreasing=TRUE)[1:20]
+df <- as.data.frame(colData(dds)[, -ncol(colData(dds))])
+pheatmap(assay(ntd)[select,], annotation_col=df)
+```
+
+
+## Principle component analysis plot
+
+```{r}
+vsd <- vst(dds, blind=FALSE)
+p = plotPCA(vsd, intgroup=c(opt$X_D)) + 
+  scale_color_discrete(name = 'Group') +
+  theme_classic()
+ggplotly(p)
+```
+