---
title: "DiffBind: CCBR ChIP-seq pipeline"
output:
    html_document:
        toc: true
        toc_depth: 2
params:
    csvfile: samplesheet.csv
    contrast: "group1_vs_group2"
    tool: "macs"
---

<!--
source: https://github.com/CCBR/Pipeliner/blob/86c6ccaa3d58381a0ffd696bbf9c047e4f991f9e/Results-template/Scripts/DiffBind_pipeliner.Rmd
-->

```{r, include=FALSE, warning=FALSE, message=FALSE}
dateandtime <- format(Sys.time(), "%a %b %d %Y - %X")

csvfile <- params$csvfile
contrasts <- params$contrast
peakcaller <- params$tool
```

### **Groups being compared:**
####    *`r contrasts`*
### **Peak sources:**
####    *`r peakcaller`*
### **Report generated:**
####    *`r dateandtime`*

```{r setup, echo=FALSE, warning=FALSE,message=FALSE}
knitr::opts_chunk$set(echo = TRUE)
suppressMessages(library(DT))
suppressMessages(library(DiffBind))
```

<br/>

## Read in sample sheet information and peak information
```{r samples, echo=FALSE, warning=FALSE,message=FALSE}
samples <- dba(sampleSheet = csvfile, peakFormat = 'bed')
consensus <- dba.peakset(samples, consensus = DBA_CONDITION)
print(samples)
```

<br/>

## Plot raw information about the peaks
### Correlation heatmap: Only peaks
```{r heatmap1, echo=FALSE, warning=FALSE,message=FALSE,out.width = "80%",fig.align="center"}
try(plot(samples, main = ""), silent = TRUE)
```

### PCA: Only peaks
```{r PCA1, echo=FALSE, warning=FALSE,message=FALSE,fig.height=5,fig.width=5,fig.align="center",fig.caption="PCA:\nOnlyPeaks"}
try(dba.plotPCA(samples, DBA_CONDITION), silent = TRUE)
```

### Overlapping peak counts
```{r Venn, echo=FALSE, warning=FALSE,message=FALSE,fig.align="center",fig.height=5,fig.width=5}
if (nrow(samples$samples) < 5) {
  dba.plotVenn(samples, 1:nrow(samples$samples))
} else {
  dba.plotVenn(samples, samples$masks[[3]])
  dba.plotVenn(samples, samples$masks[[4]])
  dba.plotVenn(consensus, consensus$masks$Consensus)
}
```

```{r peaksORsummits, echo=F}
if (grepl("narrow", samples$samples$Peaks[1])) {
  summits <- TRUE
  print("Narrow peak calling tool.")
  print("Differential peaks are 250bp upstream and downstream of the summits.")
} else if (grepl("broad", samples$samples$Peaks[1])) {
  summits <- FALSE
  print("Broad peak calling tool.")
  print("Differential peaks are consensus peaks.")
} else {
  summits <- FALSE
  print("Indeterminate peak calling tool.")
  print("Differential peaks are consensus peaks.")
}
```

## Read in bam file information under all peaks found in at least two samples
```{r DBcount, echo=FALSE, warning=FALSE,message=FALSE}
if (summits == TRUE) {
  DBdataCounts <- dba.count(samples, summits = 250)
} else {
  DBdataCounts <- dba.count(samples)
}
print(DBdataCounts)
```

<br/>

## Plot raw information about all analyzed peaks
### Correlation heatmap: Peaks and reads
```{r heatmap2, echo=FALSE, warning=FALSE,message=FALSE,out.width = "80%",fig.align="center"}
try(plot(DBdataCounts, main = ""), silent = TRUE)
```

### Heatmap: Average signal across each peak
```{r heatmap3, echo=FALSE, warning=FALSE,message=FALSE,out.width = "80%",fig.align="center"}
try(dba.plotHeatmap(DBdataCounts, correlations = FALSE), silent = TRUE)
```

### PCA: Peaks and reads
```{r PCA2, echo=FALSE, warning=FALSE,message=FALSE,fig.height=5,fig.width=5,fig.align="center"}
try(dba.plotPCA(DBdataCounts, DBA_CONDITION), silent = TRUE)
```

## Associate individual samples with the different contrasts
```{r contrast, echo=FALSE, warning=FALSE,message=FALSE}
DBdatacontrast <- dba.contrast(DBdataCounts, minMembers = 2, categories = DBA_CONDITION)
print(DBdatacontrast)
```

<br/>

## Call differential peaks using Deseq2
```{r analyze, echo=FALSE, warning=FALSE,message=FALSE}
DBAnalysisDeseq2 <- dba.analyze(DBdatacontrast, method = DBA_DESEQ2)
```

```{r, echo=FALSE, warning=FALSE,message=FALSE}
DBReportDeseq2 <- dba.report(DBAnalysisDeseq2, method = DBA_DESEQ2)
```

### PCA
```{r PCA3, echo=FALSE, warning=FALSE,message=FALSE,fig.height=5,fig.width=5,fig.align="center"}
try(dba.plotPCA(DBAnalysisDeseq2, contrast = 1, method = DBA_DESEQ2), silent = TRUE)
```


### MA plot
```{r MA, echo=FALSE, warning=FALSE,message=FALSE,fig.width=10,fig.height=4,fig.align="center"}

try(dba.plotMA(DBAnalysisDeseq2, method = DBA_DESEQ2), silent = TRUE)
```

### Volcano plot
```{r Volcano1, echo=FALSE, warning=FALSE,message=FALSE,out.width = "80%",fig.align="center"}
try(dba.plotVolcano(DBAnalysisDeseq2, method = DBA_DESEQ2), silent = TRUE)
```

### Boxplots
```{r BoxPlot, echo=FALSE, warning=FALSE,message=FALSE,fig.width=10,fig.height=4,fig.align="center"}
par(mfcol = c(1, 2))
if (length(DBReportDeseq2) > 0) {
  try(dba.plotBox(DBAnalysisDeseq2, method = DBA_DESEQ2), silent = TRUE)
} else {
  plot(0, type = "n", axes = FALSE, ann = FALSE)
}
```

## Differentially bound peaks
```{r Deseq2Report, echo=FALSE, warning=FALSE,message=FALSE}
outfile <- paste0(contrasts, "-", peakcaller, "_Diffbind_Deseq2.txt")
write.table(DBReportDeseq2, outfile, quote = F, sep = "\t", row.names = F)
DT::datatable(data.frame(DBReportDeseq2), rownames = F)
```

## R tool version information
```{r Info, echo=FALSE, message=FALSE, warning=FALSE}
sessionInfo()
```
