## 03_DE.R
## Analysis of Sequencing Data
## Vojtech Barton, 19/11/2024

###########################
### Prepare environment ###
###########################
# install packages
if (!require("BiocManager", quietly = TRUE))
  install.packages("BiocManager")

BiocManager::install("DESeq2")
BiocManager::install("apeglm")
BiocManager::install("DEGreport")

install.packages("pheatmap")
install.packages(RColorBrewer)

# load packages
library(DESeq2)
library(pheatmap)
library(RColorBrewer)
library(DEGreport)

##################
### Input data ###
##################

# specify input data
counts <- 'input_full/counts.txt'
design <- 'input_full/design.txt'

# load input data into R
counts <- read.csv(counts, sep = '\t', row.names = 1)
design <- read.csv(design, sep = '\t', row.names = 1)


##################
### Run DESeq2 ###
##################

# run DESeq2
dds <- DESeqDataSetFromMatrix(countData = counts,
                              colData = design,
                              design = ~ condition)
dds

# Pre-filtering
# While it is not necessary to pre-filter low count genes 
# before running the DESeq2 functions, there are two reasons 
# which make pre-filtering useful: by removing rows in which 
# there are very few reads, we reduce the memory size of the dds 
# data object, and we increase the speed of the transformation 
# and testing functions within DESeq2. Here we perform a minimal 
# pre-filtering to keep only rows that have at least 10 reads total.
keep <- rowSums(counts(dds)) >= 10
dds <- dds[keep,]

# Genes differentially expressed between conditions
dds <- DESeq(dds)

res <- results(dds)
res

# Log fold change shrinkage for visualization and ranking
# Shrinkage of effect size (LFC estimates) is useful for 
# visualization and ranking of genes. To shrink the LFC, 
# we pass the dds object to the function lfcShrink. Below 
# we specify to use the apeglm method for effect size shrinkage.
resultsNames(dds)

resLFC <- lfcShrink(dds, coef="condition_untreat_vs_treat", type="apeglm")
resLFC

# Ordering by smallest p-value
resOrdered <- res[order(res$pvalue),]
resOrdered

# summarizing results
summary(res)


######################
### Visualizations ###
######################

# MA-plot
plotMA(res, ylim=c(-2,2))
plotMA(resLFC, ylim=c(-2,2))

# Plot counts
plotCounts(dds, gene=which.min(res$padj), intgroup="condition")

# Heatmap
vsd <- vst(dds)
select <- order(rowMeans(counts(dds,normalized=TRUE)),
                decreasing=TRUE)[1:20]
df <- as.data.frame(colData(dds)[,c("condition","sequencing")])
pheatmap(assay(vsd)[select,], cluster_rows=FALSE, show_rownames=FALSE,
         cluster_cols=FALSE, annotation_col=df)

# Heatmap of sample-to-sample distances
sampleDists <- dist(t(assay(vsd)))
sampleDistMatrix <- as.matrix(sampleDists)
rownames(sampleDistMatrix) <- paste(vsd$condition, vsd$type, sep="-")
colnames(sampleDistMatrix) <- NULL
colors <- colorRampPalette( rev(brewer.pal(9, "Blues")) )(255)
pheatmap(sampleDistMatrix,
         clustering_distance_rows=sampleDists,
         clustering_distance_cols=sampleDists,
         col=colors)

# PCA
plotPCA(vsd, intgroup=c("condition", "sequencing"))

# Gene plots
degPlot(dds = dds, res = res, n = 6, xs = "condition")
degPlotWide(counts = dds, genes = row.names(res)[1:6], group = "condition")

# Volcano plot
degVolcano(data.frame(res[,c("log2FoldChange","padj")]))


##
resreport <- degResults(dds = dds, name = "test", org = NULL,
                        do_go = FALSE, group = "condition", xs = "condition",
                        path_results = "report")
