---
title: 'DESeq2 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 = opt$echo,
  error = TRUE
)
```


# User input

```{r 'user input'}
df = data.frame(name = names(opt)[-1],
                value = unlist(opt))
datatable(df, rownames = FALSE)
```


# Count Matrix

Display the first 100 rows of count data matrix.

```{r 'count matrix'}
count_data = read.table(opt$count_data)
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)
```

# 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')
```