---
title: "Step 6: Tidy up Count Matrix"
subtitle: "Conduct prefiltering & identify outliers"
author: "Sarah Tanja"
date: 2024-11-18
date-modified: today
format: gfm
prefer-html: true
toc: true
toc-depth: 3
link-external-icon: true
link-external-newwindow: true
reference-location: margin
citation-location: margin
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(
  echo = TRUE,         # Display code chunks
  eval = TRUE,        # Evaluate code chunks
  warning = FALSE,     # Hide warnings
  message = FALSE,     # Hide messages
  comment = ""         # Prevents appending '##' to beginning of lines in code output
)
```

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)

# Load packages

```{r}
library(tidyverse)
library(kableExtra)
```

# Import metadata

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

## Format metadata

```{r}

metadata$embryonic_stage <- gsub(" ", "", metadata$development_stage) # Remove space in development_stage condition column and rename to embryonic_stage
metadata$leachate <- metadata$pvc_leachate_concentration_mg.L # Rename to leachate to simplify

# Make a new column that represents the 12 distinct treatment groups (embryonic_phase x pvc_leachate_level) in the experiment
metadata <- metadata %>% 
  mutate(treatment_group = paste(embryonic_stage, pvc_leachate_level, sep="_")) %>% 
  relocate(treatment_group, .after = collection_date)

# Make a new column that represents the parent crosses
metadata <- metadata %>% 
  mutate(parents = sub("^(.*?)[CLMH].*", "\\1", sample_name)) %>% 
  relocate(parents, .after = sample_name)

# Make a new column that simplifies the treatment group name
metadata <- metadata %>% 
  mutate(group = sub(".*?([CLMH].*)", "\\1", sample_name)) %>% 
  relocate(group, .after = parents)

# Reorder factors
metadata$embryonic_stage <- as.factor(metadata$embryonic_stage)
metadata$embryonic_stage <- fct_relevel(metadata$embryonic_stage, "cleavage", "prawnchip", "earlygastrula")
 
metadata$pvc_leachate_level <- as.factor(metadata$pvc_leachate_level)
metadata$pvc_leachate_level <- fct_relevel(metadata$pvc_leachate_level, "control", "low", "mid", "high")

metadata$hpf <- as.factor(metadata$hours_post_fertilization)
metadata$hpf <- fct_relevel(metadata$hpf, "4", "9", "14")
```

## Filter metadata

```{r}
metadata <- metadata %>% dplyr::select(sample_no, sample_name, parents, group, treatment_group, hpf, embryonic_stage, leachate, hours_post_fertilization, pvc_leachate_level)

metadata_filt <- metadata %>% 
  filter(!sample_name %in% c("131415L4", "789C4"))

str(metadata)
```

## Write `metadata_filt.csv`
```{r}
#output metadata for other analyses 
write_csv(metadata_filt, "../metadata/metadata_filt.csv")
```

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

# 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 probabilistically 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

# Remove outliers 
Two samples had poor hisat alignment and should be removed from the dataset: 131415L4 and 789C4

```{r}
gcm_filt <- gcm %>% 
  dplyr::select(-`131415L4`, -`789C4`)

nrow(gcm)
nrow(gcm_filt)
```

# Prefiltering

[Prefiltering](https://bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#pre-filtering)

Pre-filtering our dataset reduces the memory size of the dataframe, increases the speed of the transformation and testing functions, and improves quality of statistical analysis by removing low-coverage counts. Removed counts could represent outliers in the data and removing these improves sensitivity of statistical tests.

-   **`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)`).

### All

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

nrow(gcm)
nrow(gcm_filt)

write_csv(gcm_filt, file = "../output/06_tidyup/gcm_filt.csv")
```

::: callout-note
We reduced our gene count matrix from 54,384 to 22,634 genes by selecting genes present (in counts of at least 10!) in at least 8% of the samples
:::

We are going to run DGE analysis on each hpf separately (as well as together) so let's setup gene count matrices for each of 4, 9, and 14 hours post fertilization

### 4hpf

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

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

nrow(gcm_4)
nrow(gcm_filt_4)

save(gcm_filt_4, file = "../output/06_tidyup/gcm_filt_4.RData")
write_csv(gcm_filt_4, file = "../output/06_tidyup/gcm_filt_4.csv")
```

::: callout-note
54384 genes were filtered down to 14,191 genes that have at least 10 copies in at least 25% of the 4hpf samples (25% is representing one of the four treatment groups)
:::

### 9hpf

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

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

nrow(gcm_9)
nrow(gcm_filt_9)

save(gcm_filt_9, file = "../output/06_tidyup/gcm_filt_9.RData")
write_csv(gcm_filt_9, file = "../output/06_tidyup/gcm_filt_9.csv")
```

::: callout-note
54384 genes were filtered down to 17,368 genes that have at least 10 copies in at least 25% of the 9hpf samples (25% is representing one of the four treatment groups)
:::

### 14hpf

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

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

nrow(gcm_14)
nrow(gcm_filt_14)

save(gcm_filt_14, file = "../output/06_tidyup/gcm_filt_14.RData")
write_csv(gcm_filt_14, file = "../output/06_tidyup/gcm_filt_14.csv")
```

::: callout-note
54384 genes were filtered down to 21,447 genes that have at least 10 copies in at least 25% of the 14hpf samples (25% is representing one of the four treatment groups)
:::

::: callout-important
More gene transcripts are present as embryonic development progresses!

4hpf : 14,191

9hpf : 17,368

14hpf : 21,447
:::



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

### All
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_filt$sample_name)
gene_matrix <- (colnames(gcm_filt))
meta_check <- data.frame(meta, gene_matrix)

all(meta == gene_matrix)
```

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

```{r}
nrow(metadata_filt)
```

### 4hpf

```{r}

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

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

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

There should be 18 cleavage (4hpf) samples

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

```{r}
levels(metadata_4$leachate)
```
save metadata_4 
```{r}
#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_filt %>% 
  filter(str_detect(sample_name, "[CLMH]9"))

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

```{r}
length(metadata_9$sample_name)
length(colnames(gcm_filt_9))
```
There should be 21 prawn chip (9hpf) samples

save metadata_9 
```{r}
#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_filt %>% 
  filter(str_detect(sample_name, "[CLMH]14"))

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

```{r}
length(metadata_14$sample_name)
length(colnames(gcm_filt_14))
```
There should be 22 early gastrula (14hpf) samples

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