---
title: "05-deseq2"
output: html_document
date: "2024-02-02"
---
Rmd to run  `DESeq2` comparisons between libraries. Summer 2021 SSWD Challenge Experiments coelomocyte RNAseq libraries aligned to the gene list FASTA from the published _Pycnopodia helianthoides_ genome. 

Load packages needed to run code:
```{r}
library(DESeq2)
library(tidyverse)
library(pheatmap)
library(RColorBrewer)
library(data.table)
library(ggplot2)
```

# load rounded count matrix of all 32 coelomocyte RNAseq libraries with genome gene list FASTA alignment: 
```{r}
countmatrix <- read.delim("../data/kallisto_count_matrix_rounded.tab", header = TRUE, sep = '\t')
head(countmatrix)
```

# Subset the libraries from "Experiment B" - stars that were part of the final experiment that were injected with their inoculate (either unfiltered raw tissue homogenate from a wasting adult _Pycnopodia helianthoides_, or a heat-killed raw tissue homogenate from a wasting adult _Pycnopodia helianthoides_). Stars H03 and H05 had two time points as controls, so I picked one of each so that it's a balanced 8 vs 8 comparison.      

| library_ID | star_ID | previous_exposure           | treatment         | sample_date | experiment_day |
|------------|---------|-----------------------------|-------------------|-------------|----------------|
| PSC.56     | E06     | 0.45 live inoculate 9/23/21 | control (10/5/21) | 10/14/21    | 9              |
| PSC.52     | E11     | 0.45 live inoculate 9/23/21 | control (10/5/21) | 10/14/21    | 9              |
| PSC.54     | E18     | 0.45 live inoculate 9/23/21 | control (10/5/21) | 10/14/21    | 9              |
| PSC.61     | E09     | 0.45 live inoculate 9/23/21 | control (10/5/21) | 10/15/21    | 10             |
| PSC.64     | H05     | NA                          | control (10/5/21) | 10/15/21    | 10             |
| PSC.73     | H03     | NA                          | control (10/5/21) | 10/16/21    | 11             |
| PSC.76     | E10     | 0.45 live inoculate 9/23/21 | control (10/5/21) | 10/17/21    | 12             |
| PSC.81     | E07     | 0.45 live inoculate 9/23/21 | control (10/5/21) | 10/18/21    | 13             |
| PSC.59     | H01     | NA                          | exposed (10/5/21) | 10/14/21    | 9              |
| PSC.57     | H06     | NA                          | exposed (10/5/21) | 10/14/21    | 9              |
| PSC.69     | H04     | NA                          | exposed (10/5/21) | 10/15/21    | 10             |
| PSC.67     | H09     | NA                          | exposed (10/5/21) | 10/15/21    | 10             |
| PSC.71     | H18     | NA                          | exposed (10/5/21) | 10/15/21    | 10             |
| PSC.75     | H08     | NA                          | exposed (10/5/21) | 10/16/21    | 11             |
| PSC.78     | H10     | NA                          | exposed (10/5/21) | 10/17/21    | 12             |
| PSC.83     | H07     | NA                          | exposed (10/5/21) | 10/18/21    | 13             |

```{r}
expB <- select(countmatrix, "PSC.56", "PSC.52", "PSC.54", "PSC.61", "PSC.64", "PSC.73", "PSC.76", "PSC.81", "PSC.59", "PSC.57", "PSC.69", "PSC.67", "PSC.71", "PSC.75", "PSC.78", "PSC.83")
head(expB)
```
26581 rows, 16 columns 

Make a data frame for the comparison: 
```{r}
colData <- data.frame(condition=factor(c("control", "control", "control", "control", "control", "control", "control", "control", "exposed", "exposed", "exposed", "exposed", "exposed", "exposed", "exposed", "exposed")),
                      type=factor(rep("paired-end",16)))
rownames(colData) <- colnames(expB)
dds <- DESeqDataSetFromMatrix(countData = expB,
                              colData = colData,
                              design = ~ condition)
```


```{r}
dds <- DESeq(dds)
res <- results(dds)
res <- res[order(rownames(res)), ]
```

From [Bioconductor `DESeq2` Vignette](http://bioconductor.org/packages/release/bioc/vignettes/DESeq2/inst/doc/DESeq2.html):     
```{r}
vsd <- vst(dds, blind=FALSE)
rld <- rlog(dds, blind=FALSE)
head(assay(vsd), 3)
```

```{r}
plot <- plotPCA(vsd, intgroup=c("condition"))
nudge <- position_nudge(y = 4)
plot + geom_text(aes(label = name), position = nudge)
```

```{r}
# Count number of hits with adjusted p-value less then 0.05
dim(res[!is.na(res$padj) & res$padj <= 0.05, ])
```
7834 DEGs

```{r}
control_v_exposed_fig <- res
# The main plot
plot(control_v_exposed_fig$baseMean, control_v_exposed_fig$log2FoldChange, pch=20, cex=0.45, ylim=c(-15, 15), log="x", col="darkgray",
     main="Control vs Exposed (pval </= 0.05)",
     xlab="mean of normalized counts",
     ylab="Log2 Fold Change")
# Getting the significant points and plotting them again so they're a different color
control_v_exposed_fig.sig <- res[!is.na(res$padj) & res$padj <= 0.05, ]
points(control_v_exposed_fig.sig$baseMean, control_v_exposed_fig.sig$log2FoldChange, pch=20, cex=0.45, col="red")
# 2 FC lines
abline(h=c(-1,1), col="blue")
```

Write out list of those 7834 DEGs:
Wrote out 2024-02-02. 
```{r}
#write.table(control_v_exposed_fig.sig, "../analyses/05-deseq2/expB_DEGlist_control_v_exposed.tab", sep = "\t", row.names = T, quote = FALSE, col.names = TRUE)
```


# CLEAR ENVIRONMENT: Next - Experiment A Stars at Arm Drop vs Control/Healthy 
## load rounded count matrix of all 32 coelomocyte RNAseq libraries with genome gene list FASTA alignment: 
```{r}
countmatrix <- read.delim("../data/kallisto_count_matrix_rounded.tab", header = TRUE, sep = '\t')
head(countmatrix)
```

## Subset libraries from Experiment A: 
| library_ID | star_ID | previous_exposure | treatment      | sample_date | experiment_day |
|------------|---------|-------------------|----------------|-------------|----------------|
| PSC.19     | H16     | tank_exposure     | na             | 9/23/21     | na             |
| PSC.23     | H15     | tank_exposure     | control        | 9/26/21     | 3              |
| PSC.24     | E13     | tank_exposure     | 0.45live       | 9/27/21     | 4              |
| PSC.34     | H12     | tank_exposure     | 0.45live       | 10/2/21     | 9              |
| PSC.35     | E05     | na                | unfilteredlive | 10/2/21     | 9              |
| PSC.36     | E04     | na                | unfilteredlive | 10/2/21     | 9              |
| PSC.37     | E01     | na                | unfilteredlive | 10/2/21     | 9              |
| PSC.38     | E12     | na                | 0.45live       | 10/2/21     | 9              |
| PSC.39     | H05     | na                | control        | 10/2/21     | 9              |
| PSC.40     | H04     | na                | control        | 10/2/21     | 9              |
| PSC.42     | H03     | na                | control        | 10/3/21     | 9              |
| PSC.43     | E02     | na                | unfilteredlive | 10/3/21     | 10             |
| PSC.48     | E14     | na                | 0.45live       | 10/7/21     | 12             |
| PSC.49     | E16     | na                | 0.45live       | 10/7/21     | 12             |

```{r}
expA <- select(countmatrix, "PSC.19", "PSC.23", "PSC.24", "PSC.34", "PSC.35", "PSC.36", "PSC.37", "PSC.38", "PSC.39", "PSC.40", "PSC.42", "PSC.43", "PSC.48", "PSC.49")
head(expA)
```

Make a data frame for the comparison:    
"Exposed" will mean stars that have dropped arms, and "Control" will mean stars that are healthy, regardless of their treatment group as specified in above sample table. 
```{r}
colData <- data.frame(condition=factor(c("exposed", "exposed", "exposed", "exposed", "exposed", "exposed", "exposed", "control", "control", "control", "control", "exposed", "control", "control")),
                      type=factor(rep("paired-end",14)))
rownames(colData) <- colnames(expA)
dds <- DESeqDataSetFromMatrix(countData = expA,
                              colData = colData,
                              design = ~ condition)
```

```{r}
dds <- DESeq(dds)
res <- results(dds)
res <- res[order(rownames(res)), ]
```

From [Bioconductor `DESeq2` Vignette](http://bioconductor.org/packages/release/bioc/vignettes/DESeq2/inst/doc/DESeq2.html):     
```{r}
vsd <- vst(dds, blind=FALSE)
rld <- rlog(dds, blind=FALSE)
head(assay(vsd), 3)
```


```{r}
plot <- plotPCA(vsd, intgroup=c("condition"))
nudge <- position_nudge(y = 4)
plot + geom_text(aes(label = name), position = nudge)
```

```{r}
# Count number of hits with adjusted p-value less then 0.05
dim(res[!is.na(res$padj) & res$padj <= 0.05, ])
```
4093 DEGs

```{r}
control_v_exposed_fig <- res
# The main plot
plot(control_v_exposed_fig$baseMean, control_v_exposed_fig$log2FoldChange, pch=20, cex=0.45, ylim=c(-15, 15), log="x", col="darkgray",
     main="Control vs Exposed (pval </= 0.05)",
     xlab="mean of normalized counts",
     ylab="Log2 Fold Change")
# Getting the significant points and plotting them again so they're a different color
control_v_exposed_fig.sig <- res[!is.na(res$padj) & res$padj <= 0.05, ]
points(control_v_exposed_fig.sig$baseMean, control_v_exposed_fig.sig$log2FoldChange, pch=20, cex=0.45, col="red")
# 2 FC lines
abline(h=c(-1,1), col="blue")
```

Write out list of those 4093 DEGs:
Wrote out 2024-02-09. 
```{r}
#write.table(control_v_exposed_fig.sig, "../analyses/05-deseq2/expA_DEGlist_armdrop-exposed_v_healthy-control.tab", sep = "\t", row.names = T, quote = FALSE, col.names = TRUE)
```

# Experiment A: Pairwise Comparisons Across 3 treatment groups:
## Experiment A: Unfiltered Live vs. 0.45live inoculate (uvf)
```{r}
uvf <- select(countmatrix, "PSC.24", "PSC.34", "PSC.35", "PSC.36", "PSC.37", "PSC.38", "PSC.43", "PSC.48", "PSC.49")
head(uvf)
```
9 libraries 

Make a data frame for the comparison: 
0.45live will be designated as "microbef" for microbial fraction, and unfiltered will be "exposed"
```{r}
colDatauvf <- data.frame(condition=factor(c("microbef", "microbef", "exposed", "exposed", "exposed", "microbef", "exposed", "microbef", "microbef")),
                      type=factor(rep("paired-end", 9)))
rownames(colDatauvf) <- colnames(uvf)
ddsuvf <- DESeqDataSetFromMatrix(countData = uvf,
                              colData = colDatauvf,
                              design = ~ condition)
```


```{r}
ddsuvf <- DESeq(ddsuvf)
resuvf <- results(ddsuvf)
resuvf <- resuvf[order(rownames(resuvf)), ]
```
From [Bioconductor `DESeq2` Vignette](http://bioconductor.org/packages/release/bioc/vignettes/DESeq2/inst/doc/DESeq2.html):     
```{r}
vsduvf <- vst(ddsuvf, blind=FALSE)
rlduvf <- rlog(ddsuvf, blind=FALSE)
head(assay(vsduvf), 3)
```
```{r}
plotuvf <- plotPCA(vsduvf, intgroup=c("condition"))
nudge <- position_nudge(y = 4)
plotuvf + geom_text(aes(label = name), position = nudge)
```
Run PERMANOVA:
```{r}
vsduvf.transpose <- t(assay(vsduvf))
head(vsduvf.transpose)
```

vsd.transpose is currently a matrix, but adonis2 needs it to be a data frame. Make it into a dataframe:
```{r}
vsduvf.transpose.df=as.data.frame(vsduvf.transpose)
head(vsduvf.transpose.df)
```

make rownames into a column called LibraryID:
```{r}
library(tibble)
vsduvf.transpose.df <- tibble::rownames_to_column(vsduvf.transpose.df, "LibraryID")
head(vsduvf.transpose.df)
```
Add columns to dataframe to explain the libraries:
```{r}
vsduvf.transpose.df$condition <- c("microbef", "microbef", "exposed", "exposed", "exposed", "microbef", "exposed", "microbef", "microbef")
head(vsduvf.transpose.df)
```

run a permanova to get significance:
```{r}
permanova.uvf <- adonis2(scale(vsduvf.transpose.df[c(2:26582)]) %>%
                       replace(is.na(.), 0) ~ condition, data = vsduvf.transpose.df, method = "eu")  
permanova.uvf #Significant influence of all factors
``` 
0.367 --> NO SIG DIFFERENCE

```{r}
# Count number of hits with adjusted p-value less then 0.05
dim(resuvf[!is.na(resuvf$padj) & resuvf$padj <= 0.05, ])
```
no degs

```{r}
unfilt_v_microbef_fig <- resuvf
# The main plot
plot(unfilt_v_microbef_fig$baseMean, unfilt_v_microbef_fig$log2FoldChange, pch=20, cex=0.45, ylim=c(-15, 15), log="x", col="darkgray",
     main="Unfiltered vs Microbial Fraction (pval </= 0.05)",
     xlab="mean of normalized counts",
     ylab="Log2 Fold Change")
# Getting the significant points and plotting them again so they're a different color
unfilt_v_microbef_fig.sig <- resuvf[!is.na(resuvf$padj) & resuvf$padj <= 0.05, ]
points(unfilt_v_microbef_fig.sig$baseMean, unfilt_v_microbef_fig.sig$log2FoldChange, pch=20, cex=0.45, col="red")
# 2 FC lines
abline(h=c(-1,1), col="blue")
```


## Experiment A: Unfiltered Live vs. Control (uvc)
```{r}
uvc <- select(countmatrix, "PSC.23", "PSC.35", "PSC.36", "PSC.37", "PSC.39", "PSC.40", "PSC.42", "PSC.43")
head(uvc)
```

Make a data frame for the comparison: 
"exposed" is unfiltered homogenate, and "control" is heat-kille
```{r}
colDatauvc <- data.frame(condition=factor(c("control", "exposed", "exposed", "exposed", "control", "control", "control", "exposed")),
                      type=factor(rep("paired-end", 8)))
rownames(colDatauvc) <- colnames(uvc)
ddsuvc <- DESeqDataSetFromMatrix(countData = uvc,
                              colData = colDatauvc,
                              design = ~ condition)
```

```{r}
ddsuvc <- DESeq(ddsuvc)
resuvc <- results(ddsuvc)
resuvc <- resuvc[order(rownames(resuvc)), ]
```
From [Bioconductor `DESeq2` Vignette](http://bioconductor.org/packages/release/bioc/vignettes/DESeq2/inst/doc/DESeq2.html):     
```{r}
vsduvc <- vst(ddsuvc, blind=FALSE)
rlduvc <- rlog(ddsuvc, blind=FALSE)
head(assay(vsduvc), 3)
```

```{r}
plotuvc <- plotPCA(vsduvc, intgroup=c("condition"))
nudge <- position_nudge(y = 4)
plotuvc + geom_text(aes(label = name), position = nudge)
```

Run PERMANOVA:
```{r}
vsduvc.transpose <- t(assay(vsduvc))
head(vsduvc.transpose)
```

vsd.transpose is currently a matrix, but adonis2 needs it to be a data frame. Make it into a dataframe:
```{r}
vsduvc.transpose.df=as.data.frame(vsduvc.transpose)
head(vsduvc.transpose.df)
```

make rownames into a column called LibraryID:
```{r}
library(tibble)
vsduvc.transpose.df <- tibble::rownames_to_column(vsduvc.transpose.df, "LibraryID")
head(vsduvc.transpose.df)
```
Add columns to dataframe to explain the libraries:
```{r}
vsduvc.transpose.df$condition <- c("control", "exposed", "exposed", "exposed", "control", "control", "control", "exposed")
head(vsduvc.transpose.df)
```

run a permanova to get significance:
```{r}
permanova.uvc <- adonis2(scale(vsduvc.transpose.df[c(2:26582)]) %>%
                       replace(is.na(.), 0) ~ condition, data = vsduvc.transpose.df, method = "eu")  
permanova.uvc #Significant influence of all factors
``` 
0.222 --> no significant difference 

```{r}
# Count number of hits with adjusted p-value less then 0.05
dim(resuvc[!is.na(resuvc$padj) & resuvc$padj <= 0.05, ])
```
36 degs

```{r}
unfilt_v_control_fig <- resuvc
# The main plot
plot(unfilt_v_control_fig$baseMean, unfilt_v_control_fig$log2FoldChange, pch=20, cex=0.45, ylim=c(-15, 15), log="x", col="darkgray",
     main="Unfiltered vs Control (pval </= 0.05)",
     xlab="mean of normalized counts",
     ylab="Log2 Fold Change")
# Getting the significant points and plotting them again so they're a different color
unfilt_v_control_fig.sig <- resuvc[!is.na(resuvc$padj) & resuvc$padj <= 0.05, ]
points(unfilt_v_control_fig.sig$baseMean, unfilt_v_control_fig.sig$log2FoldChange, pch=20, cex=0.45, col="red")
# 2 FC line
abline(h=c(-1,1), col="blue")
```

Write out list of those 7834 DEGs:
Wrote out 2024-02-08. 
```{r}
#write.table(unfilt_v_control_fig.sig, "../analyses/05-deseq2/expA_DEGlist_control_v_unfiltered.tab", sep = "\t", row.names = T, quote = FALSE, col.names = TRUE)
```

## Experiment A: 0.45live vs Control (fvc)
```{r}
fvc <- select(countmatrix, "PSC.23", "PSC.24", "PSC.34", "PSC.38", "PSC.39", "PSC.40", "PSC.42", "PSC.48", "PSC.49")
head(fvc)
```

Make a data frame for the comparison: 
"exposed" is 0.45live homogenate, and "control" is heat-kille
```{r}
colDatafvc <- data.frame(condition=factor(c("control", "exposed", "exposed", "exposed", "control", "control", "control", "exposed", "exposed")),
                      type=factor(rep("paired-end", 9)))
rownames(colDatafvc) <- colnames(fvc)
ddsfvc <- DESeqDataSetFromMatrix(countData = fvc,
                              colData = colDatafvc,
                              design = ~ condition)
```


```{r}
ddsfvc <- DESeq(ddsfvc)
resfvc <- results(ddsfvc)
resfvc <- resfvc[order(rownames(resfvc)), ]
```
From [Bioconductor `DESeq2` Vignette](http://bioconductor.org/packages/release/bioc/vignettes/DESeq2/inst/doc/DESeq2.html):     
```{r}
vsdfvc <- vst(ddsfvc, blind=FALSE)
rldfvc <- rlog(ddsfvc, blind=FALSE)
head(assay(vsdfvc), 3)
```

```{r}
plotfvc <- plotPCA(vsdfvc, intgroup=c("condition"))
nudge <- position_nudge(y = 4)
plotfvc + geom_text(aes(label = name), position = nudge)
```
Run PERMANOVA:
```{r}
vsdfvc.transpose <- t(assay(vsdfvc))
head(vsdfvc.transpose)
```

vsd.transpose is currently a matrix, but adonis2 needs it to be a data frame. Make it into a dataframe:
```{r}
vsdfvc.transpose.df=as.data.frame(vsdfvc.transpose)
head(vsdfvc.transpose.df)
```

make rownames into a column called LibraryID:
```{r}
library(tibble)
vsdfvc.transpose.df <- tibble::rownames_to_column(vsdfvc.transpose.df, "LibraryID")
head(vsdfvc.transpose.df)
```
Add columns to dataframe to explain the libraries:
```{r}
vsdfvc.transpose.df$condition <- c("control", "exposed", "exposed", "exposed", "control", "control", "control", "exposed", "exposed")
head(vsdfvc.transpose.df)
```

run a permanova to get significance:
```{r}
permanova.fvc <- adonis2(scale(vsdfvc.transpose.df[c(2:26582)]) %>%
                       replace(is.na(.), 0) ~ condition, data = vsdfvc.transpose.df, method = "eu")  
permanova.fvc #Significant influence of all factors
``` 
0.846 --> no significant difference 

```{r}
# Count number of hits with adjusted p-value less then 0.05
dim(resfvc[!is.na(resfvc$padj) & resfvc$padj <= 0.05, ])
```

4 DEGs


```{r}
microbial_v_control_fig <- resfvc
# The main plot
plot(microbial_v_control_fig$baseMean, microbial_v_control_fig$log2FoldChange, pch=20, cex=0.45, ylim=c(-15, 15), log="x", col="darkgray",
     main="Microbial Fraction Vs Control (pval </= 0.05)",
     xlab="mean of normalized counts",
     ylab="Log2 Fold Change")
# Getting the significant points and plotting them again so they're a different color
microbial_v_control_fig.sig <- resfvc[!is.na(resfvc$padj) & resfvc$padj <= 0.05, ]
points(microbial_v_control_fig.sig$baseMean, microbial_v_control_fig.sig$log2FoldChange, pch=20, cex=0.45, col="red")
# 2 FC line
abline(h=c(-1,1), col="blue")
```

Write out list of those 4 DEGs:
Wrote out 2024-02-08. 
```{r}
#write.table(microbial_v_control_fig.sig, "../analyses/05-deseq2/expA_DEGlist_microbial_V_control.tab", sep = "\t", row.names = T, quote = FALSE, col.names = TRUE)
```