# Control vs High
Sarah Tanja
2025-05-08

- [Goals](#goals)
- [Load libraries](#load-libraries)
- [Filter count matrix](#filter-count-matrix)
- [Setup metadata](#setup-metadata)
- [Make DESeq Object](#make-deseq-object)
- [Visualizations](#visualizations)
  - [Count transformations](#count-transformations)
  - [PCA](#pca)
  - [Heatmap](#heatmap)
- [Run DESeq2](#run-deseq2)

# Goals

This quick analysis pairs down the complexity of multiple embryonic
stages and multiple levels of leachate to a direct comparison of our
control group to our highest exposure group. Across all stages of
embryonic development, is there any impact of our highest pvc leachate
treatment on gene expression?

# Load libraries

``` r
library(tidyverse)
library(DESeq2)
library(vsn)
library(pheatmap)
library(RColorBrewer)
```

# Filter count matrix

We want to look at the `control` group (embryos developed in 0.22 um
filtered seawater) vs. the `high` group across all embryonic stages. To
do this let’s start by filtering the count matrix to only the samples of
interest.

``` r
load("../output/07_exploration/gcm_filt_or.RData")
```

The genes are the rows. The samples are the columns.

remove low & mid

``` r
gcm_controlvshigh <- gcm_filt_or %>% 
  dplyr::select(-matches("L|M"))
```

# Setup metadata

> 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.

``` r
metadata_or <- read_csv("../metadata/metadata_or.csv")
```

The samples are rows

remove low & mid

``` r
metadata_controlvshigh <- metadata_or %>% 
  filter(grepl("C|H", sample_name))
```

confirm metadata sample_name matches gene count matrix columns. 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_controlvshigh$sample_name)
gene_matrix <- (colnames(gcm_controlvshigh))
meta_check <- data.frame(meta, gene_matrix)

all(meta == gene_matrix)
```

    [1] TRUE

# Make DESeq Object

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.

``` r
#Set DESeq2 design for 1 factor with 2 levels
dds_cvh <- DESeqDataSetFromMatrix(countData = gcm_controlvshigh,
                              colData = metadata_controlvshigh, 
                              design = ~ pvc_leachate_level)
```

# Visualizations

> 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.

## Count transformations

``` r
# This is now the variance stabilized transformed DESeq Data
vsd_cvh <- vst(dds_cvh, blind=FALSE)

# This is the r-log transformed DESeq Data
rld_cvh <- rlog(dds_cvh, blind=FALSE)
```

``` r
dds_cvh <- estimateSizeFactors(dds_cvh)

df <- bind_rows(
  as_data_frame(log2(counts(dds_cvh, normalized=TRUE)[, 1:2]+1)) %>%
         mutate(transformation = "log2(x + 1)"),
  as_data_frame(assay(vsd_cvh)[, 1:2]) %>% mutate(transformation = "vst"),
  as_data_frame(assay(rld_cvh)[, 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("Control vs. High across development Count Transformation Comparisons")
```

![](08.1_deg_control_vs_high_files/figure-commonmark/unnamed-chunk-9-1.png)

``` r
meanSdPlot(assay(vsd_cvh))
```

![](08.1_deg_control_vs_high_files/figure-commonmark/unnamed-chunk-10-1.png)

``` r
meanSdPlot(assay(rld_cvh))
```

![](08.1_deg_control_vs_high_files/figure-commonmark/unnamed-chunk-11-1.png)

## PCA

``` r
# Perform PCA on rlog transformed count data
pca_cvh <- prcomp(t(assay(rld_cvh)), scale = TRUE)

pca_cvh_plot <- plotPCA(rld_cvh, intgroup=c("pvc_leachate_level"), returnData=TRUE)

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

# Make the plot
PCA_CvH_plot <- ggplot(pca_cvh_plot, aes(PC1, PC2, 
                              color=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() + 
                              geom_text(aes(label = name), size = 3) +  # Add labels from the 'name' column
                              theme_minimal()

# Display the plot
PCA_CvH_plot
```

![](08.1_deg_control_vs_high_files/figure-commonmark/unnamed-chunk-12-1.png)

``` r
# Save the plot as a png
ggsave("../output/08_deg/PCA_control_vs_high.png", plot = PCA_CvH_plot, width = 8, height = 6, dpi = 600)
```

## Heatmap

``` r
# 1. Select the top 50 genes by mean expression
top_genes <- order(rowMeans(counts(dds_cvh, normalized = TRUE)), decreasing = TRUE)[1:50]

# 2. Extract rlog-transformed matrix
rlog_matrix <- assay(rld_cvh)[top_genes, ]

# 3. Extract only pvc_leachate_level annotation
annotation_df <- as.data.frame(colData(dds_cvh)[, "pvc_leachate_level", drop = FALSE])

# Ensure it's a factor in correct order: control, high
annotation_df$pvc_leachate_level <- factor(annotation_df$pvc_leachate_level, levels = c("control", "high"))

# 4. Order samples by pvc_leachate_level
ordered_cols <- order(annotation_df$pvc_leachate_level)

# Reorder matrix and annotation
rlog_matrix_ordered <- rlog_matrix[, ordered_cols]
annotation_ordered <- annotation_df[ordered_cols, , drop = FALSE]

# 5. Plot heatmap
heatmap_colors <- colorRampPalette(rev(brewer.pal(n = 9, "RdBu")))(100)

pheatmap(
  rlog_matrix_ordered,
  cluster_rows = FALSE,
  cluster_cols = FALSE,
  show_rownames = FALSE,
  annotation_col = annotation_ordered,
  color = heatmap_colors,
  main = "Top 50 Most Expressed Genes (Control vs High)"
)
```

![](08.1_deg_control_vs_high_files/figure-commonmark/unnamed-chunk-13-1.png)

# Run DESeq2

``` r
dds_cvh_deseq <- DESeq(dds_cvh)
```

``` r
summary(results(dds_cvh_deseq))
```


    out of 22633 with nonzero total read count
    adjusted p-value < 0.1
    LFC > 0 (up)       : 0, 0%
    LFC < 0 (down)     : 0, 0%
    outliers [1]       : 0, 0%
    low counts [2]     : 1, 0.0044%
    (mean count < 0)
    [1] see 'cooksCutoff' argument of ?results
    [2] see 'independentFiltering' argument of ?results

> [!IMPORTANT]
>
> There are zero differentially expressed genes between the control and
> high treatment groups across all developmental stages.