---
title: "Step 7: Differential Expression Analysis"
subtitle: "Using `DESeq2`"
author: "Sarah Tanja"
date: 05/05/2025
format:
  gfm: default  # or html if you want to render in HTML
toc: true
toc-depth: 3
link-external-icon: true
link-external-newwindow: true
reference-location: margin
citation-location: margin
---

# Install packages

```{r}
if ("tidyverse" %in% rownames(installed.packages()) == 'FALSE') install.packages('tidyverse') 
if ("DESeq2" %in% rownames(installed.packages()) == 'FALSE') install.packages('DESeq2')
if ("RColorBrewer" %in% rownames(installed.packages()) == 'FALSE') install.packages('RColorBrewer')
if ("ComplexHeatmap" %in% rownames(installed.packages()) == 'FALSE') install.packages('ComplexHeatmap')
if ("pheatmap" %in% rownames(installed.packages()) == 'FALSE') install.packages('pheatmap')
if ("kableExtra" %in% rownames(installed.packages()) == 'FALSE') install.packages('kableExtra')
if ("viridis" %in% rownames(installed.packages()) == 'FALSE') install.packages('viridis')
if ("factoextra" %in% rownames(installed.packages()) == 'FALSE') install.packages('factoextra')
```

```{r}
if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

BiocManager::install("ComplexHeatmap")
```

# Load packages

```{r}
library(tidyverse)
library(DESeq2)
library(RColorBrewer)
library(viridis)
library(pheatmap)
library(kableExtra)
library(ComplexHeatmap)
```

# Import metadata

```{r}
metadata <- read.csv("../metadata/metadata.csv")
metadata_4 <- read.csv("../metadata/metadata_4.csv")
metadata_9<- read.csv("../metadata/metadata_9.csv")
metadata_14 <- read.csv("../metadata/metadata_14.csv")
  
```

# Import DESeq Data Objects

```{r}
# Load the DESeq2 object in a different script
load("../output/06_exploration/dds_i.RData")
load("../output/06_exploration/dds_4.RData")
load("../output/06_exploration/dds_9.RData")
load("../output/06_exploration/dds_14.RData")
```

# Run DESEq2

## Interaction

```{r}
dds_i_deseq <- DESeq(dds_i)
```

```{r}
summary(results(dds_i_deseq))
```

## 4hpf

```{r}
dds_4_deseq <- DESeq(dds_4_out)
```

```{r}
summary(results(dds_4_deseq), alpha = 0.05)
```

## 9hpf

```{r}
dds_9_deseq <- DESeq(dds_9)
```

```{r}
summary(results(dds_9_deseq), alpha = 0.05)
```

## 14hpf

```{r}
dds_14_deseq <- DESeq(dds_14)
```

```{r}
summary(results(dds_14_deseq), alpha = 0.05)
```

# [Likelihood Ratio Test](https://bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#likelihood-ratio-test),

The LRT is therefore useful for testing multiple terms at once, for example testing 3 or more levels of a factor at once, or all interactions between two variables. The LRT for count data is conceptually similar to an analysis of variance (ANOVA) calculation in linear regression, except that in the case of the Negative Binomial GLM, we use an analysis of deviance (ANODEV), where the deviance captures the difference in likelihood between a full and a reduced model.

The likelihood ratio test can be performed by specifying test="LRT" when using the DESeq function, and providing a reduced design formula, e.g. one in which a number of terms from design(dds) are removed. The degrees of freedom for the test is obtained from the difference between the number of parameters in the two models. For analyses using the likelihood ratio test (using [`nbinomLRT`](http://raven.fish.washington.edu:8787/help/library/DESeq2/help/nbinomLRT)), the p-values are determined solely by the difference in deviance between the full and reduced model formula. A single log2 fold change is printed in the results table for consistency with other results table outputs, however the test statistic and p-values may nevertheless involve the testing of one or more log2 fold changes. Which log2 fold change is printed in the results table can be controlled using the `name` argument, or by default this will be the estimated coefficient for the last element of `resultsNames(object)`.

```{r}
dds_i_lrt <- DESeq(dds_i, test="LRT")
```

```{r}
dds_4_lrt <- DESeq(dds_4, test="LRT")
```

# Wald

```{r}
dds_i_wald <- DESeq(dds_i, test="Wald")
```

```{r}
dds_4_wald <- DESeq(dds_4, test="Wald")
dds_9_wald <- DESeq(dds_9, test="Wald")
dds_14_wald <- DESeq(dds_14, test="Wald")
```

```{r}
save(dds_i_wald, file = "../output/07_deg/dds_i_wald.RData")
save(dds_4_wald, file = "../output/07_deg/dds_4_wald.RData")
save(dds_9_wald, file = "../output/07_deg/dds_9_wald.RData")
save(dds_14_wald, file = "../output/07_deg/dds_14_wald.RData")
```

```{r}
summary(results(dds_i_wald), alpha = 0.05)
```

```{r}
summary(results(dds_4_wald), alpha = 0.05)
```

```{r}
summary(results(dds_9_wald), alpha = 0.05)
```

```{r}
summary(results(dds_14_wald), alpha = 0.05)
```

> `contrast=c("condition","treated","untreated")`. Multiple results can be returned for analyses beyond a simple two group comparison, so `results` takes arguments `contrast` and `name` to help the user pick out the comparisons of interest for printing a results table. The use of the `contrast` argument is recommended for exact specification of the levels which should be compared and their order.

```{r}
resultsNames(dds_i_wald)
```

```{r}
?results
```

# Embryonic Phase

# Cleavage

## Prawnchip Low vs. Mid

```{r}
prawnchip_low_vs_mid <- results(dds_i_wald, contrast=list('embryonic_phaseprawnchip.pvc_leachate_levellow', 'embryonic_phaseprawnchip.pvc_leachate_levelmid'))

prawnchip_low_vs_mid_df <- as.data.frame(subset(prawnchip_low_vs_mid, padj<0.05))

prawnchip_low_vs_mid_df$contrast <- c('prawnchip_low_vs_mid')
prawnchip_low_vs_mid_df$gene <- rownames(prawnchip_low_vs_mid_df)
rownames(prawnchip_low_vs_mid_df) <- NULL

prawnchip_low_vs_mid_df
```

```{r}
summary(prawnchip_low_vs_mid, alpha = 0.05)
```

## Prawnchip Low vs. High

```{r}
prawnchip_low_vs_high <- results(dds_i_wald, contrast=list('embryonic_phaseprawnchip.pvc_leachate_levellow', 'embryonic_phaseprawnchip.pvc_leachate_levelhigh'))
prawnchip_low_vs_high <- as.data.frame(subset(prawnchip_low_vs_high, padj<0.05))
#NBvsB_33_deg <- as.data.frame(subset(NBvsB_33_deg, abs(log2FoldChange)>1))
prawnchip_low_vs_high$contrast <- c('prawnchip_low_vs_high')
prawnchip_low_vs_high$gene <- rownames(prawnchip_low_vs_high)
rownames(prawnchip_low_vs_high) <- NULL

prawnchip_low_vs_high
```

## Prawnchip Mid vs. High

```{r}
prawnchip_mid_vs_high <- results(dds_i_wald, contrast=list('embryonic_phaseprawnchip.pvc_leachate_levelmid', 'embryonic_phaseprawnchip.pvc_leachate_levelhigh'))

prawnchip_mid_vs_high <- as.data.frame(subset(prawnchip_mid_vs_high, padj<0.05))

prawnchip_mid_vs_high$contrast <- c('prawnchip_mid_vs_high')
prawnchip_mid_vs_high$gene <- rownames(prawnchip_mid_vs_high)
rownames(prawnchip_mid_vs_high) <- NULL

prawnchip_mid_vs_high
```

From this, we are going to filter by adjusted p value \<0.05, this is generated by the LRT test statistic.

```{r}
DEG_results <- results(dds_i_wald)
DEG_results <- as.data.frame(subset(DEG_results, padj<0.05))

DEG_results$gene <- rownames(DEG_results)
rownames(DEG_results) <- NULL
dim(DEG_results)
```

There are 16,504 DEGs out of 22,634 (\~73%) of total genes that have adjusted pvalues less than 0.05... so these are the genes that are significantly differentially expressed.

```{r}
colnames(DEG_results)
```

Results have n rows (in our case 16,504) of genes and 7 columns

1.  `baseMean` The first column, `baseMean`, is the average of the normalized count values, divided by the size factors, taken over all samples in the *DESeqDataSet*.
2.  `log2FoldChange` The column `log2FoldChange` is the effect size estimate. It tells us how much the gene's expression seems to have changed due to treatment comparison. This value is reported on a logarithmic scale to base 2: for example, a log2 fold change of 1.5 means that the gene's expression is increased by a multiplicative factor of $2^{1.5} = 2.82$
3.  `lfcSE` Of course, this log2FoldChange estimate has an uncertainty associated with it, which is available in the column `lfcSE`, the standard error estimate for the log2 fold change estimate
4.  `stat`
5.  `pvalue` *DESeq2* performs for each gene a *hypothesis test* to see whether evidence is sufficient to decide against the *null hypothesis* that there is zero effect of the treatment on the gene and that the observed difference between treatment and control was merely caused by experimental variability (i.e., the type of variability that you can expect between different samples in the same treatment group). As usual in statistics, the result of this test is reported as a *p* value, and it is found in the column `pvalue`. Remember that a *p* value indicates the probability that a fold change as strong as the observed one, or even stronger, would be seen under the situation described by the null hypothesis.
6.  `padj`
7.  `gene`

## [I ran a likelihood ratio test, but results() only gives me one comparison.](https://bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#i-ran-a-likelihood-ratio-test-but-results-only-gives-me-one-comparison.)

> "... How do I get the p values for all of the variables/levels that were removed in the reduced design?"
>
> This is explained in the help page for ?results in the section about likelihood ratio test p-values, but we will restate the answer here. When one performs a likelihood ratio test, the p values and the test statistic (the stat column) are values for the test that removes all of the variables which are present in the full design and not in the reduced design. This tests the null hypothesis that all the coefficients from these variables and levels of these factors are equal to zero.
>
> The likelihood ratio test p values therefore represent a test of all the variables and all the levels of factors which are among these variables. However, the results table only has space for one column of log fold change, so a single variable and a single comparison is shown (among the potentially multiple log fold changes which were tested in the likelihood ratio test). This is indicated at the top of the results table with the text, e.g., log2 fold change (MLE): condition C vs A, followed by, LRT p-value: '\~ batch + condition' vs '\~ batch'. This indicates that the p value is for the likelihood ratio test of all the variables and all the levels, while the log fold change is a single comparison from among those variables and levels. See the help page for results for more details.

## Conduct RLOG transformation for visualization.

```{r}
DEG_results_rlog <- dds_i[unique(DEG_results$gene)]

dim(DEG_results_rlog)

DEG_results_rlog <- rlog(DEG_results_rlog) 
```

# Heatmap of zscores

Set themes and metadata.

```{r}

legend_names_col = colnames(assay(DEG_results_rlog))

# Extract the annotation dataframe
df <- as.data.frame(colData(dds_i)[,c("pvc_leachate_level", "embryonic_phase")])

names(df)<- c("pvc_leachate_level", "embryonic_phase")
rownames(df) <- legend_names_col

levels(df$pvc_leachate_level) 
levels(df$embryonic_phase)
```

Plot heatmap of z-scores of differential expression.

```{r}

# Suppress heatmap messages
ht_opt$message <- FALSE

# Display heatmap (in RStudio or default plot device)
pheatmap(
  assay(DEG_results_rlog),
  cluster_rows = TRUE,
  show_rownames = FALSE,  # Do not show row names
  color = inferno(10),    # Inferno color palette
  show_colnames = FALSE,  # Do not show column names
  fontsize_row = 3,
  scale = "row",
  cluster_cols = TRUE,
  annotation_col = df,
  labels_col = legend_names_col,
  clustering_distance_rows = "euclidean",
  clustering_distance_cols = "euclidean",
  cutree_rows = 3,
  use_raster = TRUE  # Explicitly control rasterization
)

# Save heatmap as PNG
png("../output/07_deg/zscore_heatmap.png", width = 8, height = 6, units = "in", res = 600)
pheatmap(
  assay(DEG_results_rlog),
  cluster_rows = TRUE,
  show_rownames = FALSE,  # Do not show row names
  color = inferno(10),    # Inferno color palette
  show_colnames = FALSE,  # Do not show column names
  fontsize_row = 3,
  scale = "row",
  cluster_cols = TRUE,
  annotation_col = df,
  labels_col = legend_names_col,
  clustering_distance_rows = "euclidean",
  clustering_distance_cols = "euclidean",
  cutree_rows = 3,
  use_raster = TRUE  # Explicitly control rasterization
)
dev.off()
```

# PERMANOVA of DEGs

Conduct PERMANOVA by embryonic_phase and pvc_leachate_level

Export data for PERMANOVA test.

```{r}
test_DEG<-t(assay(DEG_results_rlog)) #export as matrix
test_DEG<-as.data.frame(test_DEG)
dim(test_DEG)
```

Rows are samples (63) and columns are z scores of the differentially expressed genes (16,504)

```{r}
#add category columns

test_DEG <- test_DEG %>%
  mutate(
    sample_name = rownames(.),
    embryonic_phase = metadata$embryonic_phase[match(sample_name, metadata$sample_name)],
    pvc_leachate_level = metadata$pvc_leachate_level[match(sample_name, metadata$sample_name)],
    .before = 1 # Place these columns before the first column
  ) 

dim(test_DEG)
```

We've added 3 columns to positions \[ ,1-3\], increasing column dim from 16504 to 16507

```{r}
test_DEG$sample_name
```

```{r}
test_DEG$embryonic_phase
```

```{r}
test_DEG$pvc_leachate_level
```

Build PERMANOVA model for DEG's that are different between temperatures.

The function fviz_eig() is part of the `factoextra` package in R.

About fviz_eig() It is used to visualize the eigenvalues (also known as scree plot) in multivariate analysis such as Principal Component Analysis (PCA), Correspondence Analysis (CA), and others. The eigenvalues represent the variance explained by each principal component or dimension.

```{r}
dim(test_DEG)

scaled_test_DEG <-prcomp(test_DEG[c(4:16507)], scale=TRUE, center=TRUE)
fviz_eig(scaled_test_DEG)

# scale data
vegan_DEG <- scale(test_DEG[c(4:16507)])

# PerMANOVA 
permanova_DEG<-adonis2(vegan_DEG ~ embryonic_phase, data = test_DEG, method='eu')
permanova_DEG
```

# PCA of zscores

```{r}
zscore_PCAdata <-
  plotPCA(
    DEG_results_rlog,
    intgroup = c("pvc_leachate_level", "embryonic_phase"),
    returnData = TRUE,
    ntop =
  )

percentVar <- round(100 * attr(zscore_PCAdata, "percentVar"))

# Make the plot
zscore_PCAplot <- ggplot(data = zscore_PCAdata,
                         aes(PC1, PC2, color = embryonic_phase, shape =
                               pvc_leachate_level)) +
  geom_point(size = 3, stroke = 1) +
  xlab(paste0("PC1: ", percentVar[1], "% variance")) +
  ylab(paste0("PC2: ", percentVar[2], "% variance")) +
  coord_fixed() + 
  theme_classic() + #Set background color
  theme(panel.border = element_rect(colour = "black", fill=NA), # Set border
                     #panel.grid.major = element_blank(), #Set major gridlines 
                     #panel.grid.minor = element_blank(), #Set minor gridlines
                     axis.line = element_line(colour = "black"), #Set axes color
        plot.background=element_blank(), 
        axis.title=element_text(size=14), 
        axis.text=element_text(size=12, color="black"))  + #Set the plot background
  theme(legend.position = ("right")); zscore_PCAplot


# Save the plot as a png
ggsave(
  "../output/07_deg/zscore_PCAplot.png",
  plot = zscore_PCAplot,
  width = 8,
  height = 6,
  dpi = 600
)

```

# Perform clustering

```{r}
clustering_sig_genes <- sigLRT_genes %>%
  arrange(padj) #genes in rows, stats in columns

matrix<-assay(gvst) #genes in rows, sample ID in columns

meta <- metadata_ordered %>%
  dplyr::select(sample_id, lifestage) #sample id and lifestage column

rownames(meta)<-meta$sample_id

meta<-meta%>%
  dplyr::select(lifestage) #sample ID in rows, one column of lifestage 

# Obtain transformed values for those significant genes
cluster_vst <- matrix[clustering_sig_genes$gene, ] #rows are DEGs and columns are sample ID
```

Run with minimum 15 genes per cluster to detect large clusters because we have many DEGs. This will take awhile to run.

```{r}
# Use the `degPatterns` function from the 'DEGreport' package to show gene clusters across sample groups
clusters <- degPatterns(cluster_vst, metadata = meta, time = "lifestage", col=NULL, plot=FALSE, minc = 15)

```

# Save gene lists for functional enrichment

## Significantly differentially expressed genes

DEG_results: List of DEG's and fold change/pvalues DEG_results_rlog: Count matrix of genes from DEG list

```{r}
dim(DEG_results)
head(DEG_results)

deg_results_list<-DEG_results%>%
  dplyr::select(gene, everything())
write_csv(deg_results_list, "../output/07_deg/deg_list.csv")

dim(DEG_results_rlog)
deg_counts<-t(assay(DEG_results_rlog))
deg_counts<-as.data.frame(deg_counts)
dim(deg_counts)
deg_counts$sample<-rownames(deg_counts)
deg_counts<-deg_counts%>%
  dplyr::select(sample, everything())
write_csv(deg_counts, "../output/07_deg/deg_counts.csv")
```

## All genes detected

List of all genes in the dataset (gvst) = 22,634 genes which will serve as a reference for functional enrichment.

```{r}
rld_i <- rlog(dds_i, blind=FALSE)
save(rld_i, file = "../output/06_exploration/rld_i.RData")
```

```{r}
dim(rld_i)
all_genes<-t(assay(rld_i))
all_genes<-as.data.frame(all_genes)
dim(all_genes)
all_genes$sample<-rownames(all_genes)
all_genes<-all_genes%>%
  dplyr::select(sample, everything())

write_csv(file="../output/07_deg/all_genes.csv", all_genes)
```

# Summary & Next Steps

## Lower the false discovery rate (alpha) ?

## [Model Technical Variation](https://bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#multi-factor-designs)

> **Experiments with many samples**: in experiments with many samples (e.g. 50, 100, etc.) it is highly likely that there will be technical variation affecting the observed counts. Failing to model this additional technical variation will lead to spurious results. Many methods exist that can be used to model technical variation, which can be easily included in the DESeq2 design to control for technical variation while estimating effects of interest. See the [RNA-seq workflow](http://www.bioconductor.org/help/workflows/rnaseqGene) for examples of using RUV or SVA in combination with DESeq2. For more details on why it is important to control for technical variation in large sample experiments, see the following [thread](https://twitter.com/mikelove/status/1513468597288452097), also archived [here](https://htmlpreview.github.io/?https://github.com/frederikziebell/science_tweetorials/blob/master/DESeq2_many_samples.html) by Frederik Ziebell.

## [Interactions](https://bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#interactions)

## Nested design?
