---
title: "Gill_TC_DGE_viz"
output: html_document
date: "2024-02-14"
---



### OA vs treatment control
```{r}

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

# CaTCulate DESeq object
dds15 <- DESeqDataSetFromMatrix(countData = GOA_TC_countmatrix,
                              colData = GOA_TC_treatmentinfo,
                              design = ~ treatment)

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

PCAdata <- plotPCA(vst(dds15), 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(dds15) >= 10) >= ncol(GOA_TC_countmatrix)/3
dds15_filt <- dds15[keep,]

```

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

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

# CaTCulate DESeq object
dds16 <- DESeqDataSetFromMatrix(countData = GOW_TC_countmatrix,
                              colData = GOW_TC_treatmentinfo,
                              design = ~ treatment)

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

PCAdata <- plotPCA(vst(dds16), 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(dds16) >= 10) >= ncol(GOW_TC_countmatrix)/3
dds16_filt <- dds16[keep,]

```

### DO vs treatment control
```{r}

# Filter data
GDO_TC_treatmentinfo <- treatmentinfo %>% filter(tissue == "G" | treatment == "control" | treatment == "DO") %>% filter(!(tissue == "F" | treatment == "OA" | treatment == "OW"))
                                                                                      GDO_TC_treatmentinfo <- GDO_TC_treatmentinfo %>%
  filter(!(treatment == "control" & day == "0")) 
write.table(GDO_TC_treatmentinfo, file = "output/DEG_lists/Gill/GDO_TC_treatmentinfo.csv", row.names = FALSE)                                                                                                      
GDO_TC_countmatrix <- subset(countmatrix_2, select=row.names(GDO_TC_treatmentinfo))
write.table(GDO_TC_countmatrix, file = "output/DEG_lists/Gill/GDO_TC_countmatrix.csv", row.names = FALSE)  
# CaTCulate DESeq object
dds17 <- DESeqDataSetFromMatrix(countData = GDO_TC_countmatrix,
                              colData = GDO_TC_treatmentinfo,
                              design = ~ treatment)

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

PCAdata <- plotPCA(vst(dds17), intgroup = c("treatment"), returnData=TRUE)
# -- note: fitType='parametric', but the dispersion trend was not well captured by the
#    function: y = a/x + b, and a local regression fit was automatically substituted.
#    specify fitType='local' or 'mean' to avoid this message next time.
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(dds17) >= 10) >= ncol(GDO_TC_countmatrix)/3
dds17_filt <- dds17[keep,]

```