---
title: "DESeq_Foot_LC"
output: html_document
date: "2024-06-14"
---
```{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
install.packages("spdep")
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)
```
### Load and clean up treatment info and gene count matrix 
```{r}
treatmentinfo <- read_csv("output/treatmentinfo_clean.csv")
gene_count_matrix_clean_og <- read_delim("output/gene_count_matrix_clean.csv", delim = ",")
rownames <- gene_count_matrix_clean$...1
gene_count_matrix_clean <- subset(gene_count_matrix_clean, select = -c(...1))

rownames(gene_count_matrix_clean) <- rownames
gene_count_matrix_clean <- data.matrix(gene_count_matrix_clean)
```

### Foot all treatments 
```{r}
F_LC_treatmentinfo <- treatmentinfo %>%
  filter( treatment %in% c("OA", "OW", "control", "DO") & tissue == "F") %>%
  filter(!(treatment == "control" & day == 3))

write.table(F_LC_treatmentinfo, file = "output/DEG_lists/Foot/F_LC_treatmentinfo.csv", row.names = FALSE)

F_LC_countmatrix <- subset(gene_count_matrix_clean, select=(F_LC_treatmentinfo$sample.1))
F_LC_countmatrix <- data.frame(F_LC_countmatrix)
write.table(F_LC_countmatrix, file = "output/DEG_lists/Foot/F_LC_countmatrix.csv")
#doesn't save correctly, need to figure this out

# CaTCulate DESeq object
dds29 <- DESeqDataSetFromMatrix(countData = F_LC_countmatrix,
                              colData = F_LC_treatmentinfo,
                              design = ~ treatment)

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

PCAdata <- plotPCA(vst(dds29), 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)+
  theme_classic()+
  scale_color_manual(labels = c("Laboratory control", "Hypoxia", "Ocean acidification", "Ocean warming"), values = c("blue", "green", "red", "orange"))+
  guides(color=guide_legend("Treatment"))+
  theme(axis.text = element_text(size =18, color = "black"))+
  theme(axis.title = element_text(size =20, color = "black"))+
  theme(legend.text = element_text(size =18, color = "black"))+
  theme(legend.title = element_text(size =20, color = "black"))+
  theme(aspect.ratio = 1)



```
### All tissue all treatments
```{r}
gene_count_matrix_clean <- read_delim("output/gene_count_matrix_clean.csv", delim = ",")

LC_treatmentinfo <- treatmentinfo %>%
  filter(!(treatment == "control" & day == 3))%>%
  filter(!(tissue == "X"))

LC_countmatrix <- subset(gene_count_matrix_clean, select=(LC_treatmentinfo$sample.1))
LC_countmatrix <- data.frame(LC_countmatrix)
write.table(LC_countmatrix, file = "output/DEG_lists/LC_countmatrix.csv")
#doesn't save correctly, need to figure this out

# CaTCulate DESeq object
dds30 <- DESeqDataSetFromMatrix(countData = LC_countmatrix,
                              colData = LC_treatmentinfo,
                              design = ~ treatment + tissue + tissue:treatment)

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

PCAdata <- plotPCA(vst(dds30), intgroup = c("tissue", "treatment"), returnData=TRUE)
percentVar <- round(100*attr(PCAdata, "percentVar")) #plot PCA of samples with all data
ggplot(PCAdata, aes(PC1, PC2, color=treatment, shape = tissue)) + 
  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)+
  theme_classic()+
  scale_color_manual(labels = c("Laboratory control", "Hypoxia", "Ocean acidification", "Ocean warming"), values = c("blue", "green", "red", "orange"))+
  scale_shape_manual(labels = c("Foot", "Gill"), values = c("circle", "triangle"))+
  guides(color=guide_legend("Treatment"), shape = guide_legend("Tissue"))+
  theme(axis.text = element_text(size =18, color = "black"))+
  theme(axis.title = element_text(size =20, color = "black"))+
  theme(legend.text = element_text(size =18, color = "black"))+
  theme(legend.title = element_text(size =20, color = "black"))+
  theme(aspect.ratio = 1)

```

### OA vs lab control
```{r}
# DEG in the foot after exposure to ocean acidification
FOA_LC_treatmentinfo <- read_delim("output/DEG_lists/Foot/FOA_LC_treatmentinfo.csv", delim = " ")

FOA_LC_countmatrix <- subset(gene_count_matrix, select=(FOA_LC_treatmentinfo$sample.1))
FOA_LC_countmatrix <- data.frame(FOA_LC_countmatrix)
write.table(FOA_LC_countmatrix, file = "output/DEG_lists/Foot/FOA_LC_countmatrix.csv")
#doesn't save correctly, need to figure this out

# 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(level = 0.95)+
  theme_classic()+
  scale_color_manual(labels = c("Lab control", "Ocean acidification"), values = c("blue", "red"))+
  guides(color=guide_legend("Treatment"))

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

FOW_LC_treatmentinfo <- read_delim("output/DEG_lists/Foot/FOW_LC_treatmentinfo.csv", delim = " ")
                                                                                                      
FOW_LC_countmatrix <- subset(gene_count_matrix, select=(FOW_LC_treatmentinfo$sample.1))
write.table(FOW_LC_countmatrix, file = "output/DEG_lists/Foot/FOW_LC_countmatrix.csv")
#Doesn't work still

# 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)+
  theme_classic()+
  scale_color_manual(labels = c("Lab control", "Ocean warming"), values = c("blue", "orange"))+
  guides(color=guide_legend("Treatment"))


# 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_delim("output/DEG_lists/Foot/FDO_LC_treatmentinfo.csv", delim = " ")
                                                                                                      
FDO_LC_countmatrix <- subset(gene_count_matrix_clean, select=(FDO_LC_treatmentinfo$sample.1))
write.csv(FDO_LC_countmatrix, file = "output/DEG_lists/Foot/FDO_LC_countmatrix.csv", row.names = TRUE)
#still comes out wrong

# 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)+
  theme_classic()+
  scale_color_manual(labels = c("Lab control", "Hypoxia"), values = c("blue", "green"))+
  guides(color=guide_legend("Treatment"))

# 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 <- read_delim("output/DEG_lists/Foot/FTC_LC_treatmentinfo.csv", delim = " ")
                                                                                 
FTC_LC_countmatrix <- subset(gene_count_matrix, select=(FTC_LC_treatmentinfo$sample.1))
write.table(FTC_LC_countmatrix, file = "output/DEG_lists/Foot/FTC_LC_countmatrix.csv")    
#still wrong
FTC_LC_treatmentinfo$day <- as.factor(FTC_LC_treatmentinfo$day)

# 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(level=0.95)+
  theme_classic()+
  scale_color_manual(labels = c("Lab control", "Treatment control"), values = c("blue", "purple"))+
  guides(color=guide_legend("Treatment"))

# 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,]

```