--- a/rmarkdown_deseq2_count_matrix.Rmd	Sat Dec 30 00:25:38 2017 -0500
+++ b/rmarkdown_deseq2_count_matrix.Rmd	Sat Dec 30 16:39:39 2017 -0500
@@ -16,24 +16,92 @@
 ```
 
 
-## User input
+# User input
 
 ```{r 'user input'}
 df = data.frame(name = names(opt)[-1],
                 value = unlist(opt))
-df
+datatable(df, rownames = FALSE)
 ```
 
 
-## Count Matrix
+# Count Matrix
+
+Display the first 100 rows of count data matrix.
 
 ```{r 'count matrix'}
 count_data = read.table(opt$count_data)
-head(count_data, 10)
+col_names = trimws(strsplit(opt$count_matrix_column_names, ',')[[1]])[1:ncol(count_data)]
+col_names = col_names[!is.na(col_names)]
+colnames(count_data)[1:length(col_names)] = col_names
+datatable(head(count_data, 100))
+```
+
+# Column Data
+
+```{r 'column data'}
+col_data = read.table(opt$col_data, 
+                      stringsAsFactors = FALSE, sep=',', header = TRUE, row.names = 1)
+datatable(col_data)
+```
+
+# Match sample names
+
+The goal of this step is to rearrange the rows of the column data matrix so that the samples rows in the count data matrix and the sample columns in the count data matrix are in the same order.
+
+```{r 'match sample names'}
+col_data = col_data[col_names, ]
+datatable(col_data)
 ```
 
-```{r 'ste[ 2'}
+# DESeqDataSet
+
+```{r 'DeseqDataSet'}
+dds = DESeqDataSetFromMatrix(countData = count_data,
+                             colData = col_data,
+                             design = formula(opt$design_formula))
+dds
+```
+
+Pre-filter low count genes
 
+```{r 'pre-filtering'}
+keep = rowSums(counts(dds)) >= 10
+dds = dds[keep, ]
+dds
+```
+
+# Differential expression analysis
+
+```{r 'differential expression analysis'}
+dds = DESeq(dds)
+# res = results(dds, contrast = c(opt$contrast_condition, opt$treatment, opt$control))
+res = results(dds)
+resultsNames(dds)
+if(nrow(res) > 500) {
+  cat('The result table has more than 500 rows. Only 500 rows are randomly selected to dispaly.')
+  datatable(as.data.frame(res)[sample(1:nrow(res), 500), ])
+} else {
+  datatable(as.data.frame(res))
+}
 ```
 
 
+
+```{r 'write results into csv'}
+#Write results into a CSV file.
+write.csv(res, 'differential_genes.csv')
+```
+
+# MAplot
+
+```{r}
+plotMA(res)
+```
+
+
+```{r 'save R objects'}
+# Save R objects to a file
+save(dds, opt, file = 'deseq2.RData')
+```
+