deseqReport.Rmd.in

---
title: "PiGx-RNAseq - DESeq2 Report"
author: "BIMSB Bioinformatics Platform"
date: '`r format(as.POSIXct(if ("" != Sys.getenv("SOURCE_DATE_EPOCH")) { as.numeric(Sys.getenv("SOURCE_DATE_EPOCH")) } else { Sys.time() }, origin="1970-01-01"), "%Y-%m-%d %H:%M:%S")`'
params:
  countDataFile: ''
  colDataFile: ''
  gtfFile: ''
  description: ''
  caseSampleGroups: ''
  controlSampleGroups: ''
  covariates: ''
  prefix: ''
  workdir: '.'
  organism: ''
  logo: ''
---

<style>
#logo
{
    position: relative;
}
#logo img {
    /*position: relative;*/
    top: 25px;
    /*right: 0px;*/
    left: 50px;
    position: fixed;
    width: 125px;
    }
body
{
    position: absolute;
    top: 150px;
}
</style>

<div id="logo" align="top">
```{r echo=FALSE}
knitr::include_graphics(params$logo)
```
</div>


```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, message = FALSE, warning = FALSE)
knitr::opts_knit$set(root.dir = params$workdir)
library(ggplot2)
library(ggpubr)
library(ggrepel)
library(DESeq2)
library(DT)
library(pheatmap)
library(corrplot)
library(reshape2)
library(plotly)
library(scales)
library(crosstalk)
library(gprofiler2)
library(rtracklayer)
library(SummarizedExperiment)
# set default theme for graphics
ggplot2::theme_set(ggpubr::theme_pubclean())
```

# Description

```{r printAnalysisDescription, results='asis'}

if (!is.null(params$description)) {
   cat("  \n")
   cat(params$description)
   cat("  \n")
   cat("  \n")
}

```

PiGx RNAseq performs differential expression analysis using _DESeq2_, and produces this report. The report includes tables and figures summarizing similarities and differences between comparison groups as specified in the settings file. In addition to the differential expression analysis, there are plots and statistics for quality control of the experiment in general with an emphasis on the reproducibility of the sequencing results among the biological replicates.

This report was generated with PiGx RNAseq version @VERSION@.

# Input Settings
To reproduce the here presentined findings, consider rerunning the respective code blocks. These blocks need to be run in the order in which these appear in this file. It starts with the assignments made in the settings file for this analysis.

```{r useInputSettingsFromParamsList}

countDataFile <- params$countDataFile
colDataFile <- params$colDataFile
gtfFile <- params$gtfFile
description <- params$description
caseSampleGroups <- params$caseSampleGroups
controlSampleGroups <- params$controlSampleGroups
covariates <- params$covariates
prefix <- params$prefix
workdir <- params$workdir
organism <- gsub('\\s*', '', params$organism)
```

If the parameters list happens not to be available in the environment
you are executing this code, you may prefer to set these values directly:
```{r manualSettings,results='asis',echo = FALSE}
cat(paste("countDataFile <- \"",      countDataFile,      "\"",sep=""),
    paste("colDataFile <- \"",        colDataFile,        "\"",sep=""),
    paste("gtfFile <- \"",            gtfFile,            "\"",sep=""),
    paste("caseSampleGroups <- \"",   caseSampleGroups,   "\"",sep=""),
    paste("controlSampleGroups <- \"",controlSampleGroups,"\"",sep=""),
    paste("covariates <- \"",         covariates,         "\"",sep=""),
    paste("prefix <- \"",             prefix,             "\"",sep=""),
    paste("workdir <- \"",            workdir,            "\"",sep=""),
    paste("organism <- \"",           organism,           "\"",sep=""), sep="  \n")
```


```{r printInputSettings}
#whether to do GO analysis or not
runGO <- TRUE

#if organism is not provided, it is not possible to do GO analysis
if(organism == '') {
  runGO <- FALSE
}

#create a folder to save high quality images generated by the report script
imagesDir <- file.path(workdir, paste0(prefix, '_images'))
if(! dir.exists(imagesDir)) {
  dir.create(path = imagesDir)
}

inputParameterDesc <- c('Count Data File',
                     'Experiment Data File',
                     'GTF File',
                     'Case sample groups',
                     'Control sample groups',
                     'Covariates to control for',
                     'Prefix for output files',
                     'Working directory',
                     'Analyzed organism'
                     )
inputParameterValues <- c(countDataFile,
                          colDataFile,
                          gtfFile,
                          caseSampleGroups,
                          controlSampleGroups,
                          covariates,
                          prefix,
                          workdir,
                          organism)
inputSettings <- data.frame(parameters = inputParameterDesc,
                            values = inputParameterValues,
                            stringsAsFactors = FALSE)
DT::datatable(data = inputSettings,
              extensions = 'FixedColumns',
              options = list(fixedColumns = TRUE,
                         scrollX = TRUE,
                         pageLength = 9,
                         dom = 't'))
```

```{r prepare_inputs_import_GTF}
gtfData <- rtracklayer::import.gff(con = gtfFile, format = 'gtf')
caseSamples <- gsub("^\\s+|\\s+$", '', unlist(strsplit(x = caseSampleGroups, split = ',')))
controlSamples <- gsub("^\\s+|\\s+$", '', unlist(strsplit(x = controlSampleGroups, split = ',')))
covariates <- gsub("^\\s+|\\s+$", '', unlist(strsplit(x = covariates, split = ',')))
#read colData and countData files
colData = read.table(colDataFile, header=T, row.names = 1, sep='\t', stringsAsFactors = T, check.names = FALSE)
countData = read.table(countDataFile, header=TRUE, row.names=1, sep='\t', stringsAsFactors = T, check.names = FALSE)

#subset colData and countData - only keep case and control samples
colData <- colData[colData$group %in% c(caseSamples, controlSamples),]
countData <- subset(countData, select = rownames(colData))

#split samples as case/control for deseq
colData$AnalysisGroup <- 'Control'
colData[colData$group %in% caseSamples,]$AnalysisGroup <- 'Case'
```

```{r run_deseq2}

mapIdsToNames <- function(ids, gtfData) {
  #first figure out if the given ids are transcript or gene ids
  transcripts <- gtfData[gtfData$type == 'transcript']
  df <- unique(data.frame('transcript_id' = transcripts$transcript_id,
                   'gene_id' = transcripts$gene_id,
                   'gene_name' = transcripts$gene_name, stringsAsFactors = FALSE))
  m <- apply(head(df[,1:2], 1000), 2, function(x) sum(x %in% ids))
  #then map the ids to gene names
  if(m['transcript_id'] > m['gene_id']){
    return(df[match(ids, df$transcript_id),]$gene_name)
  } else {
    return(df[match(ids, df$gene_id),]$gene_name)
  }
}

if(length(covariates) > 0){
  designFormula <- paste("~", paste(covariates, collapse = ' + '), "+ AnalysisGroup")
} else {
  designFormula <- "~ AnalysisGroup"
}

require(DESeq2, quietly=TRUE)
message("design formula:", designFormula)
dds <- DESeq2::DESeqDataSetFromMatrix(countData = countData, colData = colData, design = stats::as.formula(designFormula))
dds <- dds[ rowSums(DESeq2::counts(dds)) > 1, ]
dds <- DESeq2::DESeq(dds)
norm.counts = DESeq2::counts(dds, normalized = TRUE)

DEtable = DESeq2::results(dds, contrast = c("AnalysisGroup", 'Case', 'Control'))
DEtable <- DEtable[order(DEtable$padj),]
DE <- as.data.frame(DEtable)

DE$geneName <- mapIdsToNames(rownames(DE), gtfData)

DEnormalizedCountsFile <- file.path(workdir, paste0(prefix, '.normalized_counts.tsv'))
write.table(x = norm.counts,
            file = DEnormalizedCountsFile,
            quote = FALSE, sep = '\t')

DEresultsFile <- file.path(workdir, paste0(prefix, '.deseq_results.tsv'))
write.table(x = DE,
            file = DEresultsFile,
            quote = FALSE, sep = '\t')

```

# Differential Expression Analysis

Differential expression (DE) analysis was done using the [DESeq2](https://bioconductor.org/packages/release/bioc/html/DESeq2.html) R package. First, read counts are transformed using a _variance stabilizing transformation_, and then the expression values of each gene is compared between the control and sample groups using a negative binomial distribution as a model.

## Differential Expression Results Table (top 1000)

This is the table of top 1000 differentially expressed genes comparing cases to controls (as specified in the input settings listed above). The table is first filtered by absolute log2foldChange > 1 and sorted by adjusted P value after multiple testing correction (`padj`). The `baseMean` refers to expression level in the controls, and the log fold change column denotes the expression in the cases, as compared to the control.

The full table of DESeq2 results and normalized counts tables can be found at:

 - **DESeq2 results table**:
 `r DEresultsFile`
 - **DESeq2 normalized counts table**:
 `r DEnormalizedCountsFile`

```{r write_DEtable}
DEsubset <- DE[!is.na(DE$padj) & DE$padj < 0.05,]
DEsubset <- DEsubset[order(DEsubset$padj),]

DT::datatable(DEsubset[1:min(nrow(DEsubset), 1000),],
          extensions = c('Buttons', 'FixedColumns', 'Scroller'),
          options = list(fixedColumns = TRUE,
                         scrollY = 400,
                         scrollX = TRUE,
                         scroller = TRUE,
                         dom = 'Bfrtip',
                         buttons = c('colvis', 'copy', 'print', 'csv','excel', 'pdf'),
                         columnDefs = list(
                           list(targets = c(3,4,5), visible = FALSE)
                         )),
          filter = 'bottom'
          )
```


# Diagnostic Plots

This section holds a number of plots meant for a quick diagnostic and/or sanity check of the analysis.

## Number of reads assigned to genes

This plot shows the number of reads, in each sample, that are assigned to genes/transcripts. Outlier samples may be faulty and should be examined.

```{r plot_readcounts}

require(ggplot2,quietly=TRUE)
require(ggrepel,quietly=TRUE)
readCounts <- as.data.frame(colSums(countData))
readCounts$group <- colData[rownames(readCounts),]$group
readCounts$sample <- rownames(readCounts)
colnames(readCounts)[1] <- 'readCounts'

quantiles <- quantile(readCounts$readCounts, c(1:20)/20)[c(1,5,15,19)]

p <- ggplot(readCounts, aes(x = sample, y = readCounts)) + geom_bar(aes(fill = group), stat = 'identity') +
  geom_hline(yintercept = as.numeric(quantiles), color = 'red') +
  geom_label_repel(data = data.frame(x = 0, y = as.numeric(quantiles)), aes(x = x, y = y, label = names(quantiles))) + theme(legend.position = 'bottom') + scale_y_continuous(labels = scales::comma) +  coord_flip()
print(p)
#save image to folder
pdf(file = file.path(imagesDir, 'readcounts.pdf'))
print(p)
invisible(dev.off())
```

## p-value histogram

The P value distribution from the DE analysis. The expected shape depends on the expected difference / similarity between the controls and samples.

```{r plot_pvalhistogram}
p <- ggplot(data = DE, aes(x = pvalue)) + geom_histogram(bins = 100)
print(p)
#save image to folder
pdf(file = file.path(imagesDir, 'pvalue_histogram.pdf'))
print(p)
invisible(dev.off())
```

## MA plot

The MA plot gives an overview of the comparison between the two groups in the experiment. The log2 fold change for each gene is plotted on the y axis, against the average expression of that gene on the x axis.

```{r plot_MA}
DESeq2::plotMA(DEtable, main=paste("MA plot"))
#save image to folder
pdf(file = file.path(imagesDir, 'MA_plot.pdf'))
DESeq2::plotMA(DEtable, main=paste("MA plot"))
invisible(dev.off())
```

```{r computePCAplots}
plotGroups <- c(covariates, 'AnalysisGroup', 'group')
pcaPlots <- lapply(plotGroups, function(g) {
  pca <- stats::prcomp(t(scale(log(norm.counts+1))), center = TRUE)
  pcaSummary <- summary(pca)
  df <- merge(as.data.frame(pca$x), colData, by = 'row.names')
  ggplot(df, aes(x = PC1, y = PC2)) +
    geom_point(aes_string(color = g)) +
    geom_label_repel(aes(label = Row.names), size = 3) +
    labs(x = paste0('PC1 (',round(pcaSummary$importance[2, 'PC1'] * 100, 1),'%)'),
         y = paste0('PC2 (',round(pcaSummary$importance[2, 'PC2'] * 100, 1),'%)')) +
    theme_bw()
})
#save image to folder
pdf(file = file.path(imagesDir, 'pcaPlots.pdf'))
for(p in pcaPlots){
  print(p)
}
invisible(dev.off())
```

## PCA plots  {.tabset}

The 2-dimensional principal component analysis plot shows which samples group together when plotted in a reduced dimension. The 2D PCA reduced dimension conserves as much of the variance in the dataset as is possible for any 2D embedding of the data. It provides a useful birds-eye view of the data and an intuition as to which factors may drive the differences between samples or groups.

```{r plotPCA, results='asis', echo = FALSE}
for (i in 1:length(pcaPlots)) {
  cat("### ",plotGroups[i],"\n")
  print(pcaPlots[[i]])
  cat('\n\n')
}
```

## Correlation Plot

The pairwise correlation plot provides a more detailed view of which samples are more similar or different.

```{r plot_corr}
M <- stats::cor(norm.counts)
corrplot::corrplot(corr = M, order = 'hclust', method = 'square', type = 'lower', tl.srt = 45, addCoef.col = 'white')
#save image to folder
pdf(file = file.path(imagesDir, 'correlationPlot.pdf'))
corrplot::corrplot(corr = M, order = 'hclust', method = 'square', type = 'lower', tl.srt = 45, addCoef.col = 'white')
invisible(dev.off())
```

## Heatmaps

### Top 100 most highly variable genes

The heatmap below summarizes the experiment, and the apparent relationship between samples, based on the 100 highest variance genes. Each column is a sample, and each row is a gene. Both rows and columns are clustered using euclidean distance and complete linkage.

```{r plot_heatmap}
select <- na.omit(names(sort(apply(X = norm.counts, MARGIN = 1, FUN = var),decreasing = T))[1:100])
df <- as.data.frame(colData[,c("group","AnalysisGroup")])
pheatmap::pheatmap(norm.counts[select,],
         cluster_rows=TRUE,
         scale = 'row',
         show_rownames=FALSE,
         cluster_cols=TRUE,
         annotation_col=df,
         main = 'Heatmap of the Normalized Expression Values of \n Top 100 Genes with highest variance across samples')
#save image to folder
pheatmap::pheatmap(norm.counts[select,],
         cluster_rows=TRUE,
         scale = 'row',
         show_rownames=FALSE,
         cluster_cols=TRUE,
         annotation_col=df,
         filename = file.path(imagesDir, 'heatmap.pdf'),
         main = 'Heatmap of the Normalized Expression Values of \n Top 100 Genes with highest variance across samples')
```

# Exploratory Plots and Tables

## Summary of up/down regulated genes - volcano plot

This volcano plot summarizes the differential expression landscape in the comparison between the two groups.

```{r plot_summary_volcano}
p <- ggplot(DE, aes(x = log2FoldChange, y = -log10(pvalue))) + geom_point(aes(color = padj < 0.1))
print(p)
#save image to folder
pdf(file = file.path(imagesDir, 'volcanoPlot.pdf'))
print(p)
invisible(dev.off())
```

## Summary of up/down regulated genes - bar plots

These bar plots summarizes the number of significantly upregulated/downregulated number of genes based on
different adjusted p-value (selected adjusted p-values are 0.001, 0.01, 0.05, and 0.1 - see facet headers) and log2 fold change thresholds (on the x-axis) used to define the significance levels.

```{r plot_summary_barplots}

filterUP <- function(df, log2fc = 1, p = 0.1) {nrow(df[df$log2FoldChange >= log2fc & !is.na(df$padj) & df$padj <= p,])}
filterDOWN <- function(df, log2fc = 1, p = 0.1) {nrow(df[df$log2FoldChange < -log2fc & !is.na(df$padj) & df$padj <= p,])}

pVals <- c(0.001, 0.01, 0.05, 0.1)
fcVals <- c(0:(max(DE$log2FoldChange)+1))

summary <- do.call(rbind, lapply(pVals, function(p) {
  do.call(rbind, lapply(fcVals, function(f){
    up <- filterUP(DE, f, p)
    down <- filterDOWN(DE, f, p)
    return(data.frame("log2FoldChange" = f, "padj" = p,
                      "upRegulated" = up, "downRegulated" = down))
  }))
}))

require(reshape2,quietly=TRUE)
mdata <- reshape2::melt(summary, id.vars = c('log2FoldChange', 'padj'))

p <- ggplot(mdata, aes(x = log2FoldChange, y = value)) + geom_bar(aes(fill = variable), stat = 'identity', position = 'dodge') + facet_grid(~ padj) + theme(legend.position = 'bottom', legend.title = element_blank()) + labs(title = 'Number of differentially up/down regulated genes', subtitle = 'based on different p-value and log2foldChange cut-off values')
print(p)
#save image to folder
pdf(file = file.path(imagesDir, 'up_down_regulated_genes_summary.pdf'))
print(p)
invisible(dev.off())
```

## Interactive box plots of genes with significant differential expression

This interactive plot lets you see genes' position in the volcano plot, as well as their expression levels in the cases and the controls (in the box plot on the left side). Use the search box to find genes of interest. Notice that only top 1000 genes that have an adjusted p-value less than 0.1 and absolute log2 fold change value of greater than 1 are plotted.

```{r plot_interactive_boxplots}
require(plotly,quietly=TRUE)
require(crosstalk,quietly=TRUE)
select <- rownames(DEsubset)
if(length(select) > 1) {
  expressionLevels <- reshape2::melt(norm.counts[select,])
  colnames(expressionLevels) <- c('geneId', 'sampleName', 'expressionLevel')

  expressionLevels$group <- colData[expressionLevels$sampleName,]$group
  expressionLevels$AnalysisGroup <- colData[expressionLevels$sampleName,]$AnalysisGroup

  matchIds <- match(expressionLevels$geneId, rownames(DE))
  expressionLevels$padj <- DE[matchIds,]$padj
  expressionLevels$log2FoldChange <- DE[matchIds,]$log2FoldChange

  sd <- crosstalk::SharedData$new(expressionLevels, ~geneId)

  lineplot <- plot_ly(sd, x = ~sampleName, y = ~expressionLevel) %>%
    plotly::group_by(geneId) %>%
    plotly::add_lines(text = ~geneId, hoverinfo = "text", color = ~AnalysisGroup)

  volcanoPlot <- plot_ly(sd, x = ~log2FoldChange, y = ~-log10(padj)) %>%
      plotly::add_markers(text = ~geneId, hoverinfo = "text")

  subplot(
  plot_ly(sd, y = ~expressionLevel, color = ~AnalysisGroup) %>%
      plotly::add_boxplot(),
  volcanoPlot
  ) %>% plotly::highlight(on =  'plotly_click', off = 'plotly_doubleclick', selectize = TRUE)
} else {
  cat("Couldn't detect at least two genes satisfying the p-value and fold change thresholds\n")
}

```



```{r results='asis'}
if(runGO == FALSE) {
  cat("Warning:Skipping GO analysis because `organism` option is not set in settings.yaml file\n")
} else if (curl::has_internet() == FALSE){
  #gprofiler2 tool needs internet access to work. So, go analysis module is conditional
  runGO <- FALSE
  cat("Warning:Skipping GO analysis as there is no internet connection to query https://biit.cs.ut.ee/gprofiler/\n")
}
```

```{r goAnalysisTitle, eval = runGO, results='asis', echo=FALSE}
cat("# GO Term Enrichment Analysis\n")
cat("\n The following table list GO terms for differentially expressed genes. GO term analysis was carried out using [g:Profiler](https://cran.r-project.org/web/packages/gprofiler2/) R package. Differentially expressed genes are defined as those genes with adjusted p-value of less than 0.1.")
```

```{r goAnalysis, eval = runGO}
#filter genes for differential expression
DEgenes <- rownames(DE[!is.na(DE$padj) & DE$padj < 0.05,])

go <- NULL

#search for enriched GO terms within the GO and pathway domains
if(length(DEgenes) > 0) {
  require(gprofiler2,quietly=TRUE)
  go <- gprofiler2::gost(query = DEgenes,
                             organism = organism,
                             significant = TRUE,
                             sources = c('GO', 'KEGG', 'REAC', 'CORUM'))

}

if(!is.null(go)) {
  goResults <- go[['result']]
  #order by p-value
  goResults <- goResults[order(goResults$p_value),]

  goResults <- subset(goResults, select = c('p_value', 'term_size', 'query_size',
                                            'intersection_size', 'precision', 'recall',
                                            'term_id', 'source', 'term_name',
                                            'effective_domain_size', 'source_order'))

  #save full GO term table to disk
  goResultsFile <- file.path(workdir, paste0(prefix, '.GOterms.tsv'))
  write.table(x = goResults,
              file = goResultsFile,
              quote = FALSE, sep = '\t')

  #only display top GO terms in the HTML report.
  max <- ifelse(nrow(goResults) > 1000, 1000, nrow(goResults))

  DT::datatable(goResults[1:max,],
          extensions = c('Buttons', 'FixedColumns', 'Scroller'),
          options = list(fixedColumns = TRUE,
                         scrollY = 400,
                         scrollX = TRUE,
                         scroller = TRUE,
                         dom = 'Bfrtip',
                         buttons = c('colvis', 'copy', 'print', 'csv','excel', 'pdf')
                         ),
          filter = 'bottom'
          )
} else {
  cat("Warning: Couldn't detect any enriched functional terms\n")
}
```


# Session Information
```{r sessionInfo}
sessionInfo()
```