---
title: "Step 6: Data Checks for Model Assumptions"
subtitle: "Following MarineOmics Functional Genomics Tutorial`"
author: "Sarah Tanja"
date: 12/02/2024
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
---

Resources:

[MarineOmics Fitting Multifactorial Models of Differential Expression](https://marineomics.github.io/DGE_comparison_v2.html)

# Install packages

```{r}
# List of packages
packages <- c(
  "DESeq2", "edgeR", "tidyverse", "ape", "vegan", "GGally", 
  "arrayQualityMetrics", "rgl", "adegenet", "MASS", 
  "data.table", "plyr", "lmtest", "reshape2", "Rmisc", "lmerTest"
)

# Loop through each package and install if not already installed
for (pkg in packages) {
  if (!pkg %in% rownames(installed.packages())) {
    install.packages(pkg, dependencies = TRUE)
  }
  library(pkg, character.only = TRUE)
}

```

```{r}
if ("tidyverse" %in% rownames(installed.packages()) == 'FALSE') install.packages('tidyverse')
if ("DESeq2" %in% rownames(installed.packages()) == 'FALSE') install.packages('DESeq2')
if ("edgeR" %in% rownames(installed.packages()) == 'FALSE') install.packages('edgeR')
if ("ape" %in% rownames(installed.packages()) == 'FALSE') install.packages('ape')
if ("vegan" %in% rownames(installed.packages()) == 'FALSE') install.packages('vegan')
if ("GGally" %in% rownames(installed.packages()) == 'FALSE') install.packages('GGally')
if ("arrayQualityMetrics" %in% rownames(installed.packages()) == 'FALSE') install.packages('arrayQualityMetrics')
if ("rgl" %in% rownames(installed.packages()) == 'FALSE') install.packages('rgl')
if ("adegenet" %in% rownames(installed.packages()) == 'FALSE') install.packages('adegenet')
if ("MASS" %in% rownames(installed.packages()) == 'FALSE') install.packages('MASS')
if ("data.table" %in% rownames(installed.packages()) == 'FALSE') install.packages('data.table')
if ("plyr" %in% rownames(installed.packages()) == 'FALSE') install.packages('plyr')
if ("lmtest" %in% rownames(installed.packages()) == 'FALSE') install.packages('lmtest')
if ("reshape2" %in% rownames(installed.packages()) == 'FALSE') install.packages('reshape2')
if ("Rmisc" %in% rownames(installed.packages()) == 'FALSE') install.packages('Rmisc')
if ("lmerTest" %in% rownames(installed.packages()) == 'FALSE') install.packages('lmerTest')
```

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

BiocManager::install(c("DESeq2", "vsn", "edgeR", "arrayQualityMetrics"))
```

# Load packages

```{r}
# Load packages
 
library(DESeq2)
library(edgeR)
library(tidyverse)
library(ape)
library(vegan)
library(GGally)
library(arrayQualityMetrics)
library(rgl)
library(adegenet)
library(MASS)
library(data.table)
library(plyr)
library(lmtest)
library(reshape2)
library(Rmisc)
library(lmerTest)
```

# Import Gene Count Matrix (gcm)

> In this count matrix, each row represents a gene, each column a sequenced RNA library, and the values give the estimated counts of fragments that were assigned to the respective gene in each library by `HISAT2`

-   `check.names = FALSE` removed the 'X' from in front of sample id's in the column headers

```{r}
gcm <- read.csv("../output/05_count/gene_count_matrix.csv", row.names="gene_id", check.names = FALSE)
head(gcm)
```

The count matrix should show counts of whole numbers (integers) ... since each column shows <int> as the data type we know that the counts are already integers. Each column name represents a sample that is described by parental crosses (the first string of numbers), the pollution exposure (C, L, M, H), and the embryonic phase (4, 9, 14).

# Import metadata

```{r}
metadata <- read.csv("../metadata/metadata.csv")
head(metadata)
```

The `treatment` variable represents the combined levels of an experimental replicate across all variables (`embryonic_phase` + `pvc_leachate_level`)

# Understanding experimental treatments

We have 63 samples that span three embryonic phases (20 cleavage, 21 prawn chip, 22 early gastrula)

```{r}
ephases <- table(metadata$embryonic_phase)

kable(ephases, caption = "Number of samples binned by embryonic phase", booktabs = TRUE) %>%
  kable_styling(full_width = FALSE, position = "center")
```

And 4 pollution exposure levels (18 control, 15 low, 15 mid, 15 high)

```{r}
pexposure <- table(metadata$pvc_leachate_level)

kable(pexposure, caption = "Number of samples binned by pollution exposure levels", booktabs = TRUE) %>%
  kable_styling(full_width = FALSE, position = "center")
```

Giving us 12 distinct treatment groups

```{r}
tgroup <- table(metadata$treatment)

# Use kable to print the table in a cleaner format
kable(tgroup, caption = "Number of samples binned by treatment", format = "markdown", booktabs = TRUE) %>%
  kable_styling(full_width = FALSE, position = "center")
```

We expect there to be different genes present in the early embryonic development phases... and we hypothesize that there will be differential expression in response the pollution exposure at each of these developmental checkpoints.

Our minimum number of samples that have a unique treatment is 5...

Therefore, we will accept genes that are present in 5 of the samples because we are hypothesizing different expression by *embryonic_phase \* pvc_leachate_level*.

# Filter and visualize read counts

Visualize Data for Negative Binomial Distribution

Data must be rounded to nearest integer in order to be fit for negative binomial distribution

```{r}
head(gcm)
```

We can see that our gene count matrix is already made of whole integers

## Plot distribution of unfiltered read counts across all samples

```{r}
ggplot(data = data.frame(rowMeans(gcm)),
       aes(x = rowMeans.gcm.)) +
  geom_histogram(fill = "grey", binwidth = 5) +
  xlim(0, 500) +
  ylim(0, 5000) +
  theme_classic() +
  labs(title = "Distribution of unfiltered reads") +
  labs(y = "Density", x = "Raw read counts",
  title = "Read count distribution: untransformed, unnormalized, unfiltered")
```

As you can see in the above plot, the raw distribution of all read counts takes on a left-skewed negative binomial distribution.

> Most DE packages assume that read counts possess a negative binomial distribution. The negative binomial distribution is an extension of distributions for binary variables such as the Poisson distribution, allowing for estimations of "equidispersion" and "overdispersion", equal and greater-than-expected variation in expression attributed to biological variability . - [MarineOmics](https://marineomics.github.io/DGE_comparison_v2.html#:~:text=library(lmerTest)-,Filter%20and%20visualize%20read%20counts,-Before%20model%20fitting)

> While looking at the distribution of your raw reads is useful, these are not in fact the data you will be inputting to tests of differential expression. Let's plot the distribution of filtered reads normalized by library size, expressed as log2 counts per million reads (logCPM). These are the reads we will use in our test, and after we plot their distribution, we will conduct one more, slightly more robust, test of our data's fit to the negative binomial distribution using residuals from fitted models. - [MarineOmics](https://marineomics.github.io/DGE_comparison_v2.html#:~:text=While%20looking%20at,from%20fitted%20models.)

## Distibution of normalized, filtered read counts

```{r}
# Make a DGEList object for edgeR
y <- DGEList(counts = gcm, remove.zeros = TRUE)

# Let's remove samples with less than 0.5 cpm (this is ~10 counts in the count file) in fewer then 5/63 samples
keep <- rowSums(cpm(y) > .5) >= 5

table(keep)
```

```{r}
# Set keep.lib.sizes = F and recalculate library sizes after filtering
y <- y[keep, keep.lib.sizes = FALSE]

y <- calcNormFactors(y)

# Calculate logCPM
df_log <- cpm(y, log = TRUE, prior.count = 2)

# Plot distribution of filtered logCPM values
ggplot(data = data.frame(rowMeans(df_log)), 
       aes(x = rowMeans.df_log.) ) +
  geom_histogram(fill = "grey") +
  theme_classic() +
  labs(y = "Density", x = "Filtered read counts (logCPM)",
       title = "Distribution of normalized, filtered read counts")
```

... The raw filtered data doesn't appear to fit a negative binomial distribution!

# MDS plot visualizing experimental factors

> Before analyzing our data, it is essential that we look at the multivariate relationships between our samples based on transcriptome-wide expression levels. Below is example code and output for a principal coordinates analysis (PCOA) plot that visualizes multifactorial RNA-seq replicates according to two predictor variables across major and minor latent variables or PCOA axes.

```{r}
# Export pcoa loadings
dds.pcoa = pcoa(vegdist(t(df_log <- cpm(y, log = TRUE, prior.count = 2)),
                          method = "euclidean") / 1000)

# Create df of MDS vector loading
scores <- dds.pcoa$vectors

## Plot pcoa loadings of each sample, grouped by hpf and pvc_leachate_level

# Calculate % variation explained by each eigenvector
percent <- dds.pcoa$values$Eigenvalues
cumulative_percent_variance <- (percent / sum( percent)) * 100

# Prepare information for pcoa plot, then plot
color <- c("steelblue1", "tomato1", "goldenrod1")
par(mfrow = c(1, 1))
plot(
  scores[, 1],
  scores[, 2],
  cex = .5,
  cex.axis = 1,
  cex.lab = 1.25,
  xlab = paste("PC1, ", round(cumulative_percent_variance[1], 2), "%"),
  ylab = paste("PC2, ", round(cumulative_percent_variance[2], 2), "%")
  )

# Add visual groupings to pcoa plot
ordihull(
  scores,
  as.factor(metadata$pvc_leachate_level),
  border = NULL,
  lty = 2,
  lwd = .5,
  label = F,
  col = color,
  draw = "polygon",
  alpha = 100,
  cex = .5
  )

ordispider(scores, as.factor(metadata$treatment), label = F) # Vectors connecting samples in same pCO2 x time group

ordilabel(scores, cex = 0.5) # Label sample IDs
```

> ... we usually identify any sample that falls outside the main group of samples by a magnitude (along PC1) of greater than 3 standard deviations. Mathematically, all that you need to do is convert your PC1 values to Z-scores and then check for those \>\|3\|. In R, get these by using prcomp() and then accessing the 'x' variable of the returned object, e.g., `pca <- prcomp(t(rna.data); pca$x`... --Kevin Blighe on [biostars forum](https://www.biostars.org/p/281767/)

::: callout-note
789C4 and 131415L4 are outliers
:::

# Non-linear effects

> One of the simplest non-linear relationships that can be fitted to the expression of a transcript across an continuous variable is a second-order polynomial, otherwise known as a quadratic function, which can be expressed as yi=μ+β1x2+β2xyi=μ+β1x2+β2x where yy is the abundance of a given transcript (ii), μμ is the intercept, and yy is the continuous variable. For the parabola generated by fitting a second-order polynomial, β1β1 \> 0 opens the parabola upwards while β1β1 \< 0 opens the parabola downwards. The vertex of the parabola is controlled by β2β2 such that when β1β1 is negative, greater β2β2 values result in the vertex falling at higher values of xx.
>
> Quadratic polynomials applied to phenotypic performance curves commonly possess negative β1β1 values with positive β2β2 values: a downard-opening parabola with a positive vertex. However, quadtratic curves fitted to gene expression data can take on a variety of postiive or negative forms similar to exponential curves, saturating curves, and parabolas. For instance, the expression of a gene may peak an intermediate level of an environmental level before crashing or it may exponentially decline across that variable. Such trends may better model changes in the transcription of a gene compared to a linear model. To get started, we will fit non-linear second order polynomials before testing for whether model predictions for a given gene are significantly improved by a non-linear linear model. - [MarineOmics](https://marineomics.github.io/DGE_comparison_v2.html#Filter_and_visualize_read_counts:~:text=for%20differential%20expression.-,Non%2Dlinear%20effects,-Gene%20expression%20traits)

# Remove outliers

from metadata

```{r}
metadata_clean <- metadata %>% 
  filter(!(sample_name == "131415L4")) %>% 
  filter(!(sample_name == "789C4"))
```

from gene count matrix

```{r}
gcm_clean <- gcm %>%
  dplyr::select(-all_of(c("131415L4", "789C4")))
```

```{r}
# Make a DGEList object for edgeR
y <- DGEList(counts = gcm_clean, remove.zeros = TRUE)

# Let's remove samples with less than 0.5 cpm (this is ~10 counts in the count file) in fewer then 5/63 samples
keep <- rowSums(cpm(y) > .5) >= 5

table(keep)
```

```{r}
# Set keep.lib.sizes = F and recalculate library sizes after filtering
y <- y[keep, keep.lib.sizes = FALSE]

y <- calcNormFactors(y)

# Calculate logCPM
df_log <- cpm(y, log = TRUE, prior.count = 2)

# Plot distribution of filtered logCPM values
ggplot(data = data.frame(rowMeans(df_log)), 
       aes(x = rowMeans.df_log.) ) +
  geom_histogram(fill = "grey") +
  theme_classic() +
  labs(y = "Density", x = "Filtered read counts (logCPM)",
       title = "Distribution of normalized, filtered read counts outliers removed")
```

```{r}
# Export pcoa loadings
dds.pcoa = pcoa(vegdist(t(df_log <- cpm(y, log = TRUE, prior.count = 2)),
                          method = "euclidean") / 1000)

# Create df of MDS vector loading
scores <- dds.pcoa$vectors

## Plot pcoa loadings of each sample, grouped by hpf and pvc_leachate_level

# Calculate % variation explained by each eigenvector
percent <- dds.pcoa$values$Eigenvalues
cumulative_percent_variance <- (percent / sum( percent)) * 100

# Prepare information for pcoa plot, then plot
color <- c("steelblue1", "tomato1", "goldenrod1")
par(mfrow = c(1, 1))
plot(
  scores[, 1],
  scores[, 2],
  cex = .5,
  cex.axis = 1,
  cex.lab = 1.25,
  xlab = paste("PC1, ", round(cumulative_percent_variance[1], 2), "%"),
  ylab = paste("PC2, ", round(cumulative_percent_variance[2], 2), "%")
  )

# Add visual groupings to pcoa plot
ordihull(
  scores,
  as.factor(metadata_clean$pvc_leachate_level),
  border = NULL,
  lty = 2,
  lwd = .5,
  label = F,
  col = color,
  draw = "polygon",
  alpha = 100,
  cex = .5
  )

ordispider(scores, as.factor(metadata_clean$treatment), label = F) # Vectors connecting samples in same pCO2 x time group

ordilabel(scores, cex = 0.5) # Label sample IDs
```

```{r}
logCPM.pca <- prcomp(t (df_log))
logCPM.pca.proportionvariances <-
((logCPM.pca$sdev ^ 2) / (sum(logCPM.pca$sdev ^ 2))) * 100

## Do treatment groups fully segregate? Wrap samples by pCO2 x time, not just pCO2
# Replot using logCPM.pca
plot(
  logCPM.pca$x,
  type = "n",
  main = NA,
  xlab = paste("PC1, ", round(logCPM.pca.proportionvariances[1], 2), "%"),
  ylab = paste("PC2, ", round(logCPM.pca.proportionvariances[2], 2), "%")
  )

points(logCPM.pca$x,
       col = "black",
       pch = 16,
       cex = 1)
       colors2 <-
       c("steelblue1",
         "dodgerblue2",
         "tomato1",
         "coral",
         "goldenrod1",
         "goldenrod3")
       
       ordihull(
       logCPM.pca$x,
       metadata_clean$treatment,
       border = NULL,
       lty = 2,
       lwd = .5,
       col = colors2,
       draw = "polygon",
       alpha = 75,
       cex = .5,
       label = T
       )
```

of 54,384 genes in the original gene count matrix (all 63 samples), which we filtered down to 46,010 after removing any genes that had 0 sequences across all of our samples.

### Tidy filter

-   **`across(everything(), ~ . > 10)`**: This applies the condition `> 10` to all columns of the data frame. For each cell, it returns `TRUE` if the value is greater than 10 and `FALSE` otherwise.

-   **`rowSums()`**: This counts the number of `TRUE` values for each row. It adds up how many columns in each row have values greater than 10.

-   **`>= 0.08 * ncol(df)`**: This checks if the count of `TRUE` values in each row is greater than or equal to 8% of the total number of columns (`ncol(df)`).

```{r}
# Filter rows where at least 8% of the columns have a value greater than 10
gcm_tidyfilt <- gcm_filt %>%
  filter(rowSums(across(everything(), ~ . > 10)) >= 0.08 * ncol(gcm))

nrow(gcm_filt)
nrow(gcm_tidyfilt)

```

Reduced from 46,010 to 22,634 genes

## Pairwise

### 4hpf

```{r}
gcm_filt_4 <- gcm %>% 
  dplyr::select(matches("[CLMH]4"))
```

```{r}
nrow(gcm_filt_4)

gcm_filt_4 <- gcm_filt_4 %>% 
  mutate(total = rowSums(across(everything()))) %>% 
  filter(total !=0) %>% 
  dplyr::select(-total)

nrow(gcm_filt_4)
```

54384 total genes filtered down to 29343 genes that were non-zero across all cleavage samples

```{r}
# Filter rows where at least (5/20) 25% of the columns have a value greater than 10
gcm_tidyfilt_4 <- gcm_filt_4 %>%
  filter(rowSums(across(everything(), ~ . > 10)) >= 0.25 * ncol(gcm_filt_4))

nrow(gcm_filt_4)
nrow(gcm_tidyfilt_4)
```

29343 non-zero genes were filtered down to 14,191 genes that have at least 10 copies in at least 25% of the cleavage samples (25% is representing one of the four treatment groups)

### 9hpf

```{r}
gcm_filt_9 <- gcm %>% 
  dplyr::select(matches("[CLMH]9"))
```

```{r}
nrow(gcm_filt_9)

gcm_filt_9 <- gcm_filt_9 %>% 
  mutate(total = rowSums(across(everything()))) %>% 
  filter(total !=0) %>% 
  dplyr::select(-total)

nrow(gcm_filt_9)
```

54384 filtered down to 38057

```{r}
# Filter rows where at least (5/20) 25% of the columns have a value greater than 10
gcm_tidyfilt_9 <- gcm_filt_9 %>%
  filter(rowSums(across(everything(), ~ . > 10)) >= 0.25 * ncol(gcm_filt_9))

nrow(gcm_filt_9)
nrow(gcm_tidyfilt_9)
```

### 14hpf

```{r}
gcm_filt_14 <- gcm %>% 
  dplyr::select(matches("[CLMH]14"))
```

```{r}
nrow(gcm_filt_14)

gcm_filt_14 <- gcm_filt_14 %>% 
  mutate(total = rowSums(across(everything()))) %>% 
  filter(total !=0) %>% 
  dplyr::select(-total)

nrow(gcm_filt_14)
```

54384 filtered down to 44058

```{r}
# Filter rows where at least (5/20) 25% of the columns have a value greater than 10
gcm_tidyfilt_14 <- gcm_filt_14 %>%
  filter(rowSums(across(everything(), ~ . > 10)) >= 0.25 * ncol(gcm_filt_14))

nrow(gcm_filt_14)
nrow(gcm_tidyfilt_14)
```

# Prep metadata for DESeq

> It is absolutely critical that the columns of the count matrix and the rows of the column data (information about samples) are in the same order. DESeq2 will not make guesses as to which column of the count matrix belongs to which row of the column data, these must be provided to DESeq2 already in consistent order.

## Full

```{r}
meta <- (metadata$sample_name)
gene_matrix <- (colnames(gcm_tidyfilt))
meta_check <- data.frame(meta, gene_matrix)

all(meta == gene_matrix)
```

There should be 63 total samples (14hpf, 9hpf, & 4hpf)

```{r}
nrow(metadata)
```

```{r}
# Render as a kableExtra table
kable(meta_check, booktabs = TRUE) %>%
  kable_styling(full_width = FALSE, position = "center")
```

## Pairwise

### 4hpf

```{r}

# Filter rows where 'sample_name' matches the pattern [CLMH]4
metadata_4 <- metadata %>% 
  filter(str_detect(sample_name, "[CLMH]4"))

# Check if all matched names are equal
all(metadata_4$sample_name == colnames(gcm_tidyfilt_4))
```

```{r}
length(metadata_4$sample_name)
length(colnames(gcm_tidyfilt_4))
```

```{r}
# Create a data frame to compare the matched column names
meta_check_4 <- data.frame(
  sample_name = metadata_4$sample_name, 
  col_name = colnames(gcm_tidyfilt_4)
)

# Render as a kableExtra table
kable(meta_check_4, booktabs = TRUE) %>%
  kable_styling(full_width = FALSE, position = "center")
```

There should be 20 cleavage (4hpf) samples

```{r}
nrow(meta_check_4)
```

`pvc_leachate_level` should be a factor with 4 levels increasing from `control`, `low`, `mid`, `high`

```{r}
levels(metadata_4$pvc_leachate_level)
```

### 9hpf

```{r}

# Filter rows where 'sample_name' matches the pattern [CLMH]9
metadata_9 <- metadata %>% 
  filter(str_detect(sample_name, "[CLMH]9"))

# Check if all matched names are equal
all(metadata_9$sample_name == colnames(gcm_tidyfilt_9))
```

```{r}
length(metadata_9$sample_name)
length(colnames(gcm_tidyfilt_9))
```

```{r}
# Create a data frame to compare the matched column names
meta_check_9 <- data.frame(
  sample_name = metadata_9$sample_name, 
  col_name = colnames(gcm_tidyfilt_9)
)

# Render as a kableExtra table
kable(meta_check_9, booktabs = TRUE) %>%
  kable_styling(full_width = FALSE, position = "center")
```

There should be 21 prawn chip (9hpf) samples

```{r}
nrow(meta_check_9)
```

### 14hpf

```{r}

# Filter rows where 'sample_name' matches the pattern [CLMH]14
metadata_14 <- metadata %>% 
  filter(str_detect(sample_name, "[CLMH]14"))

# Check if all matched names are equal
all(metadata_14$sample_name == colnames(gcm_tidyfilt_14))
```

```{r}
length(metadata_14$sample_name)
length(colnames(gcm_tidyfilt_14))
```

```{r}
# Create a data frame to compare the matched column names
meta_check_14 <- data.frame(
  sample_name = metadata_14$sample_name, 
  col_name = colnames(gcm_tidyfilt_14)
)

# Render as a kableExtra table
kable(meta_check_14, booktabs = TRUE) %>%
  kable_styling(full_width = FALSE, position = "center")
```

There should be 22 early gastrula (14hpf) samples

```{r}
nrow(meta_check_14)
```

# Make [DESeq Dataset Object from Count Matrix](https://bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#count-matrix-input)

Create a DESeqDataSet design from gene count matrix and treatment conditions. Here we set the design to look at pollution exposure (C, L, M, H) to test for any differences in gene expression across the increasing exposure concentrations within an embryonic phase.

> The simplest design formula for differential expression would be `~ condition`, where `condition` is a column in `colData(dds)` that specifies which of two (or more groups) the samples belong to. For the airway experiment, we will specify `~ cell + dex` meaning that we want to test for the effect of dexamethasone (`dex`) controlling for the effect of different cell line (`cell`).
>
> For running *DESeq2* models, you can use R's formula notation to express any fixed-effects experimental design. Note that *DESeq2* uses the same formula notation as, for instance, the *lm* function of base R. If the research aim is to determine for which genes the effect of treatment is different across groups, then interaction terms can be included and tested using a design such as `~ group + treatment + group:treatment`. See the manual page for `?results` for more examples.

## Full interaction term design

`~genotype + condition + genotype:condition`

```{r}
#Set DESeq2 design for 2 factors and their interaction
dds_i <- DESeqDataSetFromMatrix(countData = gcm_tidyfilt,
                              colData = metadata,
                              design = ~ embryonic_phase +  pvc_leachate_level + embryonic_phase:pvc_leachate_level)
```

### [Factor levels](https://bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#note-on-factor-levels)

> As a speed concern with fitting very large models, note that each additional level of a factor in the design formula adds another parameter to the GLM which is fit by DESeq2. Users might consider first removing genes with very few reads, as this will speed up the fitting procedure.

Check levels in both the `pvc_leachate_level` and `embryonic_phase` factors for the DESeq Object `dds_i` is the object made with an interaction term)

```{r}
levels(dds_i$pvc_leachate_level)
levels(dds_i$embryonic_phase)
```

Here we have 2 factors:

1.  `pvc_leachate_level` (4 levels)
    1.  `control`
    2.  `low`
    3.  `mid`
    4.  `high`
2.  `embryonic_phase` (3 levels)
    1.  `cleavage`
    2.  `prawnchip`
    3.  `earlygastrula`

## Pairwise

### 4hpf

```{r}
dds_4 <- DESeqDataSetFromMatrix(countData = gcm_tidyfilt_4,
                              colData = metadata_4,
                              design = ~ pvc_leachate_level)
```

```{r}
dds_4$pvc_leachate_level <- factor(dds_4$pvc_leachate_level, levels = c("control",
                                                                        "low",
                                                                        "mid",
                                                                        "high"))

levels(dds_4$pvc_leachate_level)
```

### 9hpf

```{r}
dds_9 <- DESeqDataSetFromMatrix(countData = gcm_tidyfilt_9,
                              colData = metadata_9,
                              design = ~ pvc_leachate_level)
```

```{r}
dds_9$pvc_leachate_level <- factor(dds_9$pvc_leachate_level, levels = c("control",
                                                                        "low",
                                                                        "mid",
                                                                        "high"))

levels(dds_9$pvc_leachate_level)
```

### 14hpf

```{r}
dds_14 <- DESeqDataSetFromMatrix(countData = gcm_tidyfilt_14,
                              colData = metadata_14,
                              design = ~ pvc_leachate_level)
```

```{r}
dds_14$pvc_leachate_level <- factor(dds_14$pvc_leachate_level, levels = c("control",
                                                                        "low",
                                                                        "mid",
                                                                        "high"))

levels(dds_14$pvc_leachate_level)
```

# Count transformations

> In order to test for differential expression, we operate on raw counts and use discrete distributions... However for visualization or clustering -- it might be useful to work with transformed versions of the count data.
>
> **Which transformation to choose?** The VST is much faster to compute and is less sensitive to high count outliers than the rlog. The rlog tends to work well on small datasets (n \< 30), potentially outperforming the VST when there is a wide range of sequencing depth across samples (an order of magnitude difference). We therefore recommend the VST for medium-to-large datasets (n \> 30). You can perform both transformations and compare the `meanSdPlot` or PCA plots generated, as described below.

## Full

```{r}
vsd_i <- vst(dds_i, blind=FALSE)
rld_i <- rlog(dds_i, blind=FALSE)
```

## 4hpf

```{r}
vsd_4 <- vst(dds_4, blind=FALSE)
rld_4 <- rlog(dds_4, blind=FALSE)
```

## 9hpf

```{r}
vsd_9 <- vst(dds_9, blind=FALSE)
rld_9 <- rlog(dds_9, blind=FALSE)
```

## 14hpf

```{r}
vsd_14 <- vst(dds_14, blind=FALSE)
rld_14 <- rlog(dds_14, blind=FALSE)
```

# [Effects of transformations on the variance](https://bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#effects-of-transformations-on-the-variance)

> To show the effect of the transformation, in the figure below we plot the first sample against the second, first simply using the *log2* function (after adding 1, to avoid taking the log of zero), and then using the VST and rlog-transformed values. For the *log2* approach, we need to first estimate *size factors* to account for sequencing depth, and then specify `normalized=TRUE`. Sequencing depth correction is done automatically for the *vst* and *rlog*.

> The figures below plot the standard deviation of the transformed data, across samples, against the mean, using the shifted logarithm transformation, the regularized log transformation and the variance stabilizing transformation. The shifted logarithm has elevated standard deviation in the lower count range, and the regularized log to a lesser extent, while for the variance stabilized data the standard deviation is roughly constant along the whole dynamic range.
>
> Note that the vertical axis in such plots is the square root of the variance over all samples, so including the variance due to the experimental conditions. While a flat curve of the square root of variance over the mean may seem like the goal of such transformations, this may be unreasonable in the case of datasets with many true differences due to the experimental conditions.

### Full

```{r}
dds_i <- estimateSizeFactors(dds_i)

df <- bind_rows(
  as_data_frame(log2(counts(dds_i, normalized=TRUE)[, 1:2]+1)) %>%
         mutate(transformation = "log2(x + 1)"),
  as_data_frame(assay(vsd_i)[, 1:2]) %>% mutate(transformation = "vst"),
  as_data_frame(assay(rld_i)[, 1:2]) %>% mutate(transformation = "rlog"))
  
colnames(df)[1:2] <- c("x", "y")  

lvls <- c("log2(x + 1)", "vst", "rlog")
df$transformation <- factor(df$transformation, levels=lvls)

ggplot(df, aes(x = x, y = y)) +
  geom_hex(bins = 80) +
  coord_fixed() +
  facet_grid(. ~ transformation) +
  ggtitle("All Samples Interaction design (4, 9 & 14 hpf) Count Transformation Comparisons")
```

```{r}
meanSdPlot(assay(vsd_i))
meanSdPlot(assay(rld_i))
```

## 4hpf

```{r}
dds_4 <- estimateSizeFactors(dds_4)

df <- bind_rows(
  as_data_frame(log2(counts(dds_4, normalized=TRUE)[, 1:2]+1)) %>%
         mutate(transformation = "log2(x + 1)"),
  as_data_frame(assay(vsd_4)[, 1:2]) %>% mutate(transformation = "vst"),
  as_data_frame(assay(rld_4)[, 1:2]) %>% mutate(transformation = "rlog"))
  
colnames(df)[1:2] <- c("x", "y")  

lvls <- c("log2(x + 1)", "vst", "rlog")
df$transformation <- factor(df$transformation, levels=lvls)

ggplot(df, aes(x = x, y = y)) +
  geom_hex(bins = 80) +
  coord_fixed() +
  facet_grid(. ~ transformation) +
  ggtitle("Cleavage Samples (4hpf) Count Transformation Comparisons")

```

```{r}
meanSdPlot(assay(vsd_4))
meanSdPlot(assay(rld_4))
```

## 9hpf

```{r}
dds_9 <- estimateSizeFactors(dds_9)

df <- bind_rows(
  as_data_frame(log2(counts(dds_9, normalized=TRUE)[, 1:2]+1)) %>%
         mutate(transformation = "log2(x + 1)"),
  as_data_frame(assay(vsd_9)[, 1:2]) %>% mutate(transformation = "vst"),
  as_data_frame(assay(rld_9)[, 1:2]) %>% mutate(transformation = "rlog"))
  
colnames(df)[1:2] <- c("x", "y")  

lvls <- c("log2(x + 1)", "vst", "rlog")
df$transformation <- factor(df$transformation, levels=lvls)

ggplot(df, aes(x = x, y = y)) +
  geom_hex(bins = 80) +
  coord_fixed() +
  facet_grid(. ~ transformation) +
  ggtitle("Prawn Chip Samples (9hpf) Count Transformation Comparisons")

```

```{r}
meanSdPlot(assay(vsd_9))
meanSdPlot(assay(rld_9))
```

## 14hpf

```{r}
dds_14 <- estimateSizeFactors(dds_14)

df <- bind_rows(
  as_data_frame(log2(counts(dds_14, normalized=TRUE)[, 1:2]+1)) %>%
         mutate(transformation = "log2(x + 1)"),
  as_data_frame(assay(vsd_14)[, 1:2]) %>% mutate(transformation = "vst"),
  as_data_frame(assay(rld_14)[, 1:2]) %>% mutate(transformation = "rlog"))
  
colnames(df)[1:2] <- c("x", "y")  

lvls <- c("log2(x + 1)", "vst", "rlog")
df$transformation <- factor(df$transformation, levels=lvls)

ggplot(df, aes(x = x, y = y)) +
  geom_hex(bins = 80) +
  coord_fixed() +
  facet_grid(. ~ transformation) +
  ggtitle("Early Gastrula Samples (14hpf) Count Transformation Comparisons")

```

```{r}
meanSdPlot(assay(vsd_14))
meanSdPlot(assay(rld_14))
```

# Identify the Outlier in 4hpf

> ... we usually identify any sample that falls outside the main group of samples by a magnitude (along PC1) of greater than 3 standard deviations. Mathematically, all that you need to do is convert your PC1 values to Z-scores and then check for those \>\|3\|. In R, get these by using prcomp() and then accessing the 'x' variable of the returned object, e.g., `pca <- prcomp(t(rna.data); pca$x`... --Kevin Blighe on [biostars forum](https://www.biostars.org/p/281767/)

```{r}
# Perform PCA on raw count matrix
pca_out_4 <- prcomp(t(gcm_tidyfilt_4), scale = TRUE)
```

```{r}
pca_out_plot_4 <- plotPCA(vsd_4, intgroup=c("pvc_leachate_level", "embryonic_phase"), returnData=TRUE)

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

# Make the plot
PCAout4 <- ggplot(pca_out_plot_4, aes(PC1, PC2, 
                              color=embryonic_phase,
                              shape=pvc_leachate_level)) +
                              geom_text(aes(label = name), size = 3) +  # Add labels from the 'name' column
                              geom_point(size=3) +
                              xlab(paste0("PC1: ",percentVar[1],"% variance")) +
                              ylab(paste0("PC2: ",percentVar[2],"% variance")) + 
                              coord_fixed() + 
                              theme_minimal()

# Display the plot
PCAout4

# Save the plot as a png
ggsave("../output/06_exploration/PCAout4.png", plot = PCAout4, width = 8, height = 6, dpi = 600)
```

::: callout-important
Samples `131415L4` and `789C4`may be outliers in the cleavage phase
:::

## Remove the outliers

### From metadata

```{r}
metadata_4_out <- metadata_4 %>% 
  filter(!(sample_name == "131415L4")) %>% 
  filter(!(sample_name == "789C4"))
```

### From gene count matrix

```{r}
gcm_tidyfilt_4_out <- gcm_tidyfilt_4 %>% 
  select(-"131415L4", -"789C4")

nrow(gcm_tidyfilt_4_out)
```

14191 genes

### Refilter

```{r}
# Filter rows where at least (3/18) = .16% of the columns have a value greater than 10
gcm_tidyfilt_4_out <- gcm_tidyfilt_4_out %>%
  filter(rowSums(across(everything(), ~ . > 10)) >= 0.16 * ncol(gcm_tidyfilt_4_out))

nrow(gcm_tidyfilt_4_out)
```

14191 genes (same genes present with or without outlier samples)

### Remake DESeq Object

```{r}
dds_4_out <- DESeqDataSetFromMatrix(countData = gcm_tidyfilt_4_out,
                              colData = metadata_4_out,
                              design = ~ pvc_leachate_level)
```

### Redo count transformations

```{r}
vsd_4_out <- vst(dds_4_out, blind=FALSE)
rld_4_out <- rlog(dds_4_out, blind=FALSE)
```

```{r}
dds_4_out <- estimateSizeFactors(dds_4_out)

df <- bind_rows(
  as_data_frame(log2(counts(dds_4_out, normalized=TRUE)[, 1:2]+1)) %>%
         mutate(transformation = "log2(x + 1)"),
  as_data_frame(assay(vsd_4_out)[, 1:2]) %>% mutate(transformation = "vst"),
  as_data_frame(assay(rld_4_out)[, 1:2]) %>% mutate(transformation = "rlog"))
  
colnames(df)[1:2] <- c("x", "y")  

lvls <- c("log2(x + 1)", "vst", "rlog")
df$transformation <- factor(df$transformation, levels=lvls)

ggplot(df, aes(x = x, y = y)) +
  geom_hex(bins = 80) +
  coord_fixed() +
  facet_grid(. ~ transformation) +
  ggtitle("Cleavage Samples (4hpf) Count Transformation Comparisons Outliers Removed")

```

```{r}
meanSdPlot(assay(vsd_4_out))
title(main = "VST of 4hpf with outliers removed")
meanSdPlot(assay(rld_4_out))
```

## Remove outlier samples from Full dataset

### From metadata

```{r}
metadata_out <- metadata %>% 
  filter(!(sample_name == "131415L4")) %>% 
  filter(!(sample_name == "789C4"))
```

### From gene count matrix

```{r}
gcm_tidyfilt_out <- gcm_tidyfilt %>% 
  select(-"131415L4", -"789C4")

nrow(gcm_tidyfilt_out)
```

22634 genes

### Refilter

```{r}
# Filter rows where at least (3/60) = .05% of the columns have a value greater than 10
gcm_tidyfilt_out <- gcm_tidyfilt_out %>%
  filter(rowSums(across(everything(), ~ . > 10)) >= 0.05 * ncol(gcm_tidyfilt_out))

nrow(gcm_tidyfilt_out)
```

22634 genes (same genes present with or without outlier samples)

### Remake DESeq Object

```{r}
dds_i_out <- DESeqDataSetFromMatrix(countData = gcm_tidyfilt_out,
                              colData = metadata_out,
                              design = ~ embryonic_phase +  pvc_leachate_level + embryonic_phase:pvc_leachate_level)
```

### Redo count transformations

```{r}
vsd_i_out <- vst(dds_i_out, blind=FALSE)
rld_i_out <- rlog(dds_i_out, blind=FALSE)
```

```{r}
dds_i_out <- estimateSizeFactors(dds_i_out)

df <- bind_rows(
  as_data_frame(log2(counts(dds_i_out, normalized=TRUE)[, 1:2]+1)) %>%
         mutate(transformation = "log2(x + 1)"),
  as_data_frame(assay(vsd_i_out)[, 1:2]) %>% mutate(transformation = "vst"),
  as_data_frame(assay(rld_i_out)[, 1:2]) %>% mutate(transformation = "rlog"))
  
colnames(df)[1:2] <- c("x", "y")  

lvls <- c("log2(x + 1)", "vst", "rlog")
df$transformation <- factor(df$transformation, levels=lvls)

ggplot(df, aes(x = x, y = y)) +
  geom_hex(bins = 80) +
  coord_fixed() +
  facet_grid(. ~ transformation) +
  ggtitle("All Samples Count Transformation Comparisons Outliers Removed")

```

```{r}
# Extract PC1 values
pc1_values <- pca_out_4$x[, "PC1"]

# Convert PC1 values to Z-scores
pc1_zscores <- scale(pc1_values)

# Identify samples with |Z| > 3
outliers <- abs(pc1_zscores) > 3

# Filter for outlier samples
outlier_samples <- which(outliers)

# Output results
outlier_data <- data.frame(
  Sample = rownames(scaled_test$x)[outlier_samples],
  PC1 = pc1_values[outlier_samples],
  Z_score = pc1_zscores[outlier_samples]
)

print(outlier_data)
```

# ADD PERMANOVA of FULL "Big Picture" GCM

Build a distance matrix. To transform a gene count matrix into a distance matrix, you need to decide on an appropriate distance or dissimilarity measure. Common options include:

-   Bray-Curtis dissimilarity: Often used for ecological and compositional data.
-   Euclidean distance: Measures straight-line distance and may not be ideal for count data.
-   Manhattan distance: Sum of absolute differences.
-   Jaccard index: Based on presence/absence of genes.

Export data for PERMANOVA test.

```{r}
test<-t(assay(rld_i)) #export as matrix
test<-as.data.frame(test)
test
```

This is a matrix of rlog transformed count data for all genes across all embryonic phases, prior to DGE analysis

```{r}
#add category columns
test <- test %>% 
  mutate(
    sample_name = rownames(test),
    embryonic_phase = metadata$embryonic_phase[match(test$sample_name, metadata$sample_name)],
    pvc_leachate_level = metadata$pvc_leachate_level[match(test$sample_name, metadata$sample_name)]
        ) %>% 
  relocate(sample_name, embryonic_phase, pvc_leachate_level, .before = 1)

test
```

## Build PERMANOVA model & PCA

```{r}
# Perform PCA
scaled_test <-prcomp(test[c(4:22637)], scale=TRUE, center=TRUE)

# Calculate proportion of variance explained
eigenvalues <- scaled_test$sdev^2
variance_explained <- eigenvalues / sum(eigenvalues) * 100

# Create a data frame for plotting
eig_df <- data.frame(
  Component = seq_along(variance_explained),
  Variance = variance_explained
)

# Plot with ggplot2
ggplot(eig_df, aes(x = Component, y = Variance)) +
  geom_bar(stat = "identity", fill = "skyblue") +
  geom_text(aes(label = round(Variance, 1)), vjust = -0.5) +
  labs(
    title = "Scree Plot",
    x = "Principal Components",
    y = "Percentage of Variance Explained"
  ) +
  theme_minimal()

# from factoextra package
fviz_eig(scaled_test)
```

```{r}
test$embryonic_phase <- as.factor(test$embryonic_phase)
test$pvc_leachate_level <- as.factor(test$pvc_leachate_level)

levels(test$embryonic_phase)
levels(test$pvc_leachate_level)
```

## Run PERMANOVA

```{r}
# scale data
vegan <- scale(test[c(4:22637)])

# Test embryonic_phase alone
permanova_embryo <- adonis2(
  vegan ~ embryonic_phase,
  data = test,
  method = "euclidean"
)
print(permanova_embryo)

# Test pvc_leachate_level alone
permanova_pvc <- adonis2(
  vegan ~ pvc_leachate_level,
  data = test,
  method = "euclidean"
)
print(permanova_pvc)

# Test interaction
permanova_interaction <- adonis2(
  vegan ~ embryonic_phase:pvc_leachate_level,
  data = test,
  method = "euclidean"
)
print(permanova_interaction)


# PerMANOVA sequential
permanova <- adonis2(
  vegan ~ embryonic_phase + pvc_leachate_level + embryonic_phase:pvc_leachate_level,
  data = test,
  method = "euclidean"
)

print(permanova)
```

## NMDS

```{r}
nmds_result <- metaMDS(vegan, distance = "euclidean")
plot(nmds_result, display = "sites")
```

## Betadisper

> One measure of multivariate dispersion (variance) for a group of samples is to calculate the average distance of group members to the group centroid or spatial median (both referred to as 'centroid' from now on unless stated otherwise) in multivariate space. To test if the dispersions (variances) of one or more groups are different, the distances of group members to the group centroid are subject to ANOVA. This is a multivariate analogue of Levene's test for homogeneity of variances if the distances between group members and group centroids is the Euclidean distance.
>
> However, better measures of distance than the Euclidean distance are available for ecological data. These can be accommodated by reducing the distances produced using any dissimilarity coefficient to principal coordinates, which embeds them within a Euclidean space. The analysis then proceeds by calculating the Euclidean distances between group members and the group centroid on the basis of the principal coordinate axes rather than the original distances.

```{r}
?betadispr
```

# Exploratory Data Visualization

::: callout-caution
USE THE DESEQ OBJECT FOR EXPLORATORY DATA VIS
:::

## PCA

> A principal components analysis (PCA) plot shows the samples in the 2D plane spanned by their first two principal components. This type of plot is useful for visualizing the overall effect of experimental covariates and batch effects.

plotPCA(vsd, intgroup=c("pvc_leachate_level", "embryonic_phase"))

```{r}
vsd_pca <- plotPCA(vsd_i, intgroup=c("pvc_leachate_level", "embryonic_phase"), returnData=TRUE)

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

# Make the plot
vsdPCA <- ggplot(vsd_pca, aes(PC1, PC2, 
                              color=embryonic_phase,
                              shape=pvc_leachate_level)) +
                              geom_point(size=3) +
                              xlab(paste0("PC1: ",percentVar[1],"% variance")) +
                              ylab(paste0("PC2: ",percentVar[2],"% variance")) + 
                              coord_fixed()

# Display the plot
vsdPCA

# Save the plot as a png
ggsave("../output/06_exploration/vsdPCA.png", plot = vsdPCA, width = 8, height = 6, dpi = 600)
```

```{r}
rld_pca <- plotPCA(rld_i, intgroup=c("pvc_leachate_level", "embryonic_phase"), returnData=TRUE)

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

# Create PCA plot with labels
rldPCA <- ggplot(rld_pca, aes(PC1, PC2, 
                              color = embryonic_phase,
                              shape = pvc_leachate_level,
                              label = name)) +  # Add sample_name for labeling
  geom_point(size = 3) +  # Points with color and shape
  #geom_text_repel(size = 3) +  # Labels with repel to avoid overlap
  xlab(paste0("PC1: ", percentVar[1], "% variance")) +
  ylab(paste0("PC2: ", percentVar[2], "% variance")) +
  coord_fixed() +  # Maintain aspect ratio
  theme_minimal()  # Optional: Apply a minimal theme

# Display the plot
print(rldPCA)

# Save the plot as a png
ggsave("../output/06_exploration/rldPCA.png", plot = rldPCA, width = 8, height = 6, dpi = 600)
```

## Heatmaps

-   `normalized = TRUE` means the DESeq2 Data Set was run through `dds_i <- estimateSizeFactors(dds_i)` in the \[effects of transformations on the variance\] code chunk above

### Heatmap PVC Leachate Level

```{r}
# Select the top 20 rows based on mean expression
select <- order(rowMeans(counts(dds_i, normalized = TRUE)), decreasing = TRUE)[1:20]

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

# Ensure pvc_leachate_level is a factor with the desired order
heat$pvc_leachate_level <- factor(heat$pvc_leachate_level, levels = c("control", "low", "mid", "high"))

# Order columns by pvc_leachate_level
ordered_columns <- order(heat$pvc_leachate_level)

# Reorder the assay matrix and annotation dataframe
rld_matrix <- assay(rld_i)[select, ][, ordered_columns]
rld_heat <- heat[ordered_columns, ]

# Plot the heatmap
explore_plevel_heatmap <- pheatmap(
                            rld_matrix,
                            cluster_rows = FALSE,          # Do not cluster rows
                            show_rownames = FALSE,         # Do not show row names
                            cluster_cols = FALSE,          # Disable column clustering
                            annotation_col = rld_heat      # Add annotation
                                  )
# display the heatpmap
explore_plevel_heatmap

# Save the heatmap as a png
ggsave("../output/06_exploration/explore_plevel_heatmap.png", plot = explore_plevel_heatmap, width = 8, height = 6, dpi = 600)
```

### Heatmap Embryonic Phase

```{r}
# Select the top 20 rows based on mean expression
select <- order(rowMeans(counts(dds_i, normalized = TRUE)), decreasing = TRUE)[1:20]

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

# Ensure pvc_leachate_level is a factor with the desired order
heat$embryonic_phase <- factor(heat$embryonic_phase, levels = c("cleavage", "prawnchip", "earlygastrula"))

# Order columns by pvc_leachate_level
ordered_columns <- order(heat$embryonic_phase)

# Reorder the assay matrix and annotation dataframe
rld_matrix <- assay(rld_i)[select, ][, ordered_columns]
rld_heat <- heat[ordered_columns, ]

# Plot the heatmap
explore_ephase_heatmap <- pheatmap(
                            rld_matrix,
                            cluster_rows = FALSE,          # Do not cluster rows
                            show_rownames = FALSE,         # Do not show row names
                            cluster_cols = FALSE,          # Disable column clustering
                            annotation_col = rld_heat      # Add annotation
                                  )
# display the heatpmap
explore_ephase_heatmap

# Save the heatmap as a png
ggsave("../output/06_exploration/explore_ephase_heatmap.png", plot = explore_ephase_heatmap, width = 8, height = 6, dpi = 600)
```

# Summary & Next Steps

```{r}
save(dds_i, file = "../output/06_exploration/dds_i.RData")

#output metadata for other analyses 
write_csv(metadata, "../metadata/metadata.csv")
```

```{r}
save(dds_4, file = "../output/06_exploration/dds_4.RData")

#output metadata for other analyses 
write_csv(metadata_4, "../metadata/metadata_4.csv")
```

```{r}
save(dds_9, file = "../output/06_exploration/dds_9.RData")

#output metadata for other analyses 
write_csv(metadata_9, "../metadata/metadata_9.csv")
```

```{r}
save(dds_14, file = "../output/06_exploration/dds_14.RData")

#output metadata for other analyses 
write_csv(metadata_14, "../metadata/metadata_14.csv")
```
