---
title: "DGE_viz"
output: html_document
date: "2024-02-05"
---
```{r load libraries}
if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install(version = "3.16")

BiocManager::install("DESeq2")
library("DESeq2") # Bioconductor
library("tidyverse")
library("dplyr")
library("pheatmap")
library("RColorBrewer")
library("genefilter") # Bioconductor
library("ggplot2")
library("gplots")
library("limma") # Bioconductor
library("spdep") 
library("adegenet")
library("tximport") # Bioconductor
BiocManager::install("vsn")
BiocManager::install("tidybulk")
BiocManager::install("goseq")
BiocManager::install("affycoretools")
BiocManager::install("EnhancedVolcano")
BiocManager::install("pcaExplorer")
BiocManager::install("apeglm")
BiocManager::install("PCAtools")
BiocManager::install("DESeq2PCA")
library(pcaExplorer)
library(PCAtools)
```


## Overall comparison with treatment control
```{r}
F_treatmentinfo  <- treatmentinfo.2 %>% filter(tissue == "F") %>%
  filter(!( tissue == "G"))

F_treatmentinfo <- F_treatmentinfo %>%
  filter(!(treatment == "control" & day == "0")) 
write.table(F_treatmentinfo, file = "output/DEG_lists/Foot/F_treatmentinfo.csv", row.names = FALSE)

F_countmatrix <- subset(countmatrix_2, select=row.names(F_treatmentinfo))
write.table(F_countmatrix, file = "output/DEG_lists/Foot/F_countmatrix.csv", row.names = FALSE)

dds3 <- DESeqDataSetFromMatrix(countData = F_countmatrix,
                              colData = F_treatmentinfo,
                              design = ~ treatment)

dds3 <- DESeq(dds3)
resultsNames(dds3) # lists the coefficients
summary(dds3)

PCAdataF <- plotPCA(vst(dds3), intgroup = c("treatment"), returnData=TRUE)
percentVarF <- round(100*attr(PCAdataF, "percentVar")) #plot PCA of samples with all data
ggplot(PCAdataF, aes(PC1, PC2, color=treatment, show.legend = FALSE)) + 
  geom_point(size=4, alpha = 5/10) +
  xlab(paste0("PC1: ",percentVarF[1],"% variance")) +
  ylab(paste0("PC2: ",percentVarF[2],"% variance")) + 
  coord_fixed()+
  stat_ellipse(level=0.95, show.legend = FALSE)+
  theme_classic()+
  labs(color='Treatment')+
  scale_color_manual(labels = c("Control", "Hypoxia", "Ocean Acidification", "Ocean Warming"), values = c("blue", "purple","green", "red"))+
  theme(text=element_text(size = 20))
  
Foot
```
### Grid arrange of two overall PCA's
```{r}

install.packages("gridExtra")
library(gridExtra)
library(cowplot)

grid.arrange(Foot, Gill, nrow = 1)
plot_grid(Foot, Gill, nrow = 1, rel_heights = c(0.5, 2))

```



### OA vs lab control
```{r}
# DEG in the foot after exposure to ocean acidification

# Filter data
FOA_LC_treatmentinfo <- treatmentinfo.2 %>% filter(tissue == "F" | treatment == "control" | treatment == "OA") %>% filter(!(tissue == "G" | treatment == "OW" | treatment == "DO"))
                                                                                                      FOA_LC_treatmentinfo <- FOA_LC_treatmentinfo %>%
  filter(!(treatment == "control" & day == "3")) 
write.table(FOA_LC_treatmentinfo, file = "output/DEG_lists/Foot/FOA_LC_treatmentinfo.csv", row.names = FALSE)

FOA_LC_countmatrix <- subset(countmatrix_2, select=row.names(FOA_LC_treatmentinfo))
write.table(FOA_LC_countmatrix, file = "output/DEG_lists/Foot/FOA_LC_countmatrix.csv", row.names = FALSE)

# Calculate DESeq object
dds4 <- DESeqDataSetFromMatrix(countData = FOA_LC_countmatrix,
                              colData = FOA_LC_treatmentinfo,
                              design = ~ treatment)

dds4 <- DESeq(dds4)
resultsNames(dds4) # lists the coefficients

PCAdata <- plotPCA(vst(dds4), intgroup = c("treatment"), returnData=TRUE)
percentVar <- round(100*attr(PCAdata, "percentVar")) #plot PCA of samples with all data
ggplot(PCAdata, aes(PC1, PC2, color=treatment)) + 
  geom_point(size=4, alpha = 5/10) +
  xlab(paste0("PC1: ",percentVar[1],"% variance")) +
  ylab(paste0("PC2: ",percentVar[2],"% variance")) + 
  coord_fixed()+
  stat_ellipse()

# Filtering: keep genes that have at least 10 counts across 1/3 of the samples - https://support.bioconductor.org/p/110307/
keep <- rowSums(DESeq2::counts(dds4) >= 10) >= ncol(FOA_LC_countmatrix)/3
dds4_filt <- dds4[keep,]

```

### OW vs lab control
```{r}

# Filter data
FOW_LC_treatmentinfo <- treatmentinfo.2 %>% filter(tissue == "F" | treatment == "control" | treatment == "OW") %>% filter(!(tissue == "G" | treatment == "OA" | treatment == "DO"))
                                                                               FOW_LC_treatmentinfo <- FOW_LC_treatmentinfo %>%filter(!(treatment == "control" & day == "3")) 
write.table(FOW_LC_treatmentinfo, file = "output/DEG_lists/Foot/FOW_LC_treatmentinfo.csv", row.names = FALSE)
                                                                                                      
FOW_LC_countmatrix <- subset(countmatrix_2, select=row.names(FOW_LC_treatmentinfo))
write.table(FOW_LC_countmatrix, file = "output/DEG_lists/Foot/FOW_LC_countmatrix.csv", row.names = FALSE)

# Calculate DESeq object
dds5 <- DESeqDataSetFromMatrix(countData = FOW_LC_countmatrix,
                              colData = FOW_LC_treatmentinfo,
                              design = ~ treatment)

dds5 <- DESeq(dds5)
resultsNames(dds5) # lists the coefficients

PCAdata <- plotPCA(vst(dds5), intgroup = c("treatment"), returnData=TRUE)
percentVar <- round(100*attr(PCAdata, "percentVar")) #plot PCA of samples with all data
ggplot(PCAdata, aes(PC1, PC2, color=treatment)) + 
  geom_point(size=4, alpha = 5/10) +
  xlab(paste0("PC1: ",percentVar[1],"% variance")) +
  ylab(paste0("PC2: ",percentVar[2],"% variance")) + 
  coord_fixed()+
  stat_ellipse(level=0.95)

# Filtering: keep genes that have at least 10 counts across 1/3 of the samples - https://support.bioconductor.org/p/110307/
keep <- rowSums(DESeq2::counts(dds5) >= 10) >= ncol(FOW_LC_countmatrix)/3
dds5_filt <- dds5[keep,]

```

### DO vs lab control
```{r}
FDO_LC_treatmentinfo <- read_table("output/DEG_lists/Foot/FDO_LC_treatmentinfo.csv")


# Filter data
FDO_LC_treatmentinfo <- treatmentinfo %>% filter(tissue == "F" | treatment == "control" | treatment == "DO") %>% filter(!(tissue == "G" | tissue == "X" | treatment == "OA" | treatment == "OW"))
                                                                                                    FDO_LC_treatmentinfo <- FDO_LC_treatmentinfo %>%
  filter(!(treatment == "control" & day == "3")) 
write.table(FDO_LC_treatmentinfo, file = "output/DEG_lists/Foot/FDO_LC_treatmentinfo.csv", row.names = FALSE)
                                                                                                      
FDO_LC_countmatrix <- subset(countmatrix_2, select=row.names(FDO_LC_treatmentinfo))
write.table(FDO_LC_countmatrix, file = "output/DEG_lists/Foot/FDO_LC_countmatrix.csv", row.names = FALSE)

# Calculate DESeq object
dds6 <- DESeqDataSetFromMatrix(countData = FDO_LC_countmatrix,
                              colData = FDO_LC_treatmentinfo,
                              design = ~ treatment)

dds6 <- DESeq(dds6)
resultsNames(dds6) # lists the coefficients

PCAdata <- plotPCA(vst(dds6), intgroup = c("treatment"), returnData=TRUE)
percentVar <- round(100*attr(PCAdata, "percentVar")) #plot PCA of samples with all data
ggplot(PCAdata, aes(PC1, PC2, color=treatment)) + 
  geom_point(size=4, alpha = 5/10) +
  xlab(paste0("PC1: ",percentVar[1],"% variance")) +
  ylab(paste0("PC2: ",percentVar[2],"% variance")) + 
  coord_fixed()+
  stat_ellipse(level=0.95)

# Filtering: keep genes that have at least 10 counts across 1/3 of the samples - https://support.bioconductor.org/p/110307/
keep <- rowSums(DESeq2::counts(dds6) >= 10) >= ncol(FDO_LC_countmatrix)/3
dds6_filt <- dds6[keep,]

```


### treatment control vs lab control
```{r}
# RERUN THE DDS, accidentally overwrote it 
# Filter data
FTC_LC_treatmentinfo <- treatmentinfo.2 %>% filter(tissue == "F" | treatment == "control") %>% filter(!(tissue == "G" | treatment == "OA" | treatment == "OW" | treatment == "DO"))
 write.table(FTC_LC_treatmentinfo, file = "output/DEG_lists/Foot/FTC_LC_treatmentinfo.csv", row.names = FALSE)                                                                                                 
FTC_LC_countmatrix <- subset(countmatrix_2, select=row.names(FTC_LC_treatmentinfo))
write.table(FTC_LC_countmatrix, file = "output/DEG_lists/Foot/FTC_LC_countmatrix.csv", row.names = FALSE)    

# Calculate DESeq object
dds7 <- DESeqDataSetFromMatrix(countData = FTC_LC_countmatrix,
                              colData = FTC_LC_treatmentinfo,
                              design = ~ day)

dds7 <- DESeq(dds7)
resultsNames(dds7) # lists the coefficients

PCAdata <- plotPCA(vst(dds7), intgroup = c("day"), returnData=TRUE)
percentVar <- round(100*attr(PCAdata, "percentVar")) #plot PCA of samples with all data
ggplot(PCAdata, aes(PC1, PC2, color=day)) + 
  geom_point(size=4, alpha = 5/10) +
  xlab(paste0("PC1: ",percentVar[1],"% variance")) +
  ylab(paste0("PC2: ",percentVar[2],"% variance")) + 
  coord_fixed()+
  stat_ellipse()

# Filtering: keep genes that have at least 10 counts across 1/3 of the samples - https://support.bioconductor.org/p/110307/
keep <-  rowSums(DESeq2::counts(dds7) >= 10) >= ncol(FTC_LC_countmatrix)/3
dds7_filt <- dds7[keep,]

```