---
title: "Step 6: Data Exploration of Count Matrix"
subtitle: "Using `DESeq2`"
author: "Sarah Tanja"
date: 11/18/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:

[RNA-seq workflow](https://master.bioconductor.org/packages/release/workflows/vignettes/rnaseqGene/inst/doc/rnaseqGene.html)

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

# Install packages

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

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

BiocManager::install("DESeq2")
```

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

BiocManager::install("vsn")
```

# Load packages

```{r}
library(tidyverse)
library(DESeq2)
library(kableExtra)
library(pheatmap)
library(vsn)
library(vegan)
library(factoextra)
```

# Load metadata and filtered gene count matrices

```{r}
load("../output/06_exploration/metadata.RData")
load("../output/06_exploration/gcm_tidyfilt.RData")
load("../output/06_exploration/gcm_tidyfilt_4.RData")
load("../output/06_exploration/gcm_tidyfilt_9.RData")
load("../output/06_exploration/gcm_tidyfilt_14.RData")
```


# 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
The following code checks that all the sample names in the metadata dataframe match the column names in the gene count matrix. It should return 'TRUE'
```{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)
```
Visually check that sample names in the metadata dataframe match the column names in the gene count matrix.
```{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)
```
save metadata_4 
```{r}
save(metadata_4, file = "../metadata/metadata_4.RData")

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


### 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)
```
save metadata_9 
```{r}
save(metadata_9, file = "../metadata/metadata_9.RData")

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

### 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)
```
save metadata_14
```{r}
save(metadata_14, file = "../metadata/metadata_14.RData")

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

# 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))
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}
meanSdPlot(assay(vsd_i_out))
meanSdPlot(assay(rld_i_out))
```

# 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")
```
