---
title: "DESeq2_DEsign_Run_10cpm"
author: "Samuel Gurr"
date: "March 29, 2021"
output:
html_document:
toc: true
toc_float: true
---

# Setup: 

```{r setup, include=TRUE}
knitr::opts_chunk$set(echo = TRUE, warning = FALSE, message = FALSE, cache = TRUE)
```

### Load libraries
```{r load_libraries, include = TRUE}
# load libraries - notes show the install command needed to install (pre installed)
library(DESeq2) # note: this was previously installed with the command `BiocManager::install("DESeq2")`
library(edgeR)
library(goseq)
library(dplyr)
library(GenomicFeatures)
library(data.table)
library(calibrate)
library(affycoretools) # note: this was previously installed with the BiocManager::install("affycoretools")
library(data.table)
library(vsn)
library(tidybulk)
# Plotting
library(ggplot2)
library(cowplot)
library(pheatmap)
library(gplots)
library(RColorBrewer)
library(EnhancedVolcano)  # note: this was previously installed with the command `BiocManager::install("EnhancedVolcano")`
library(pcaExplorer) 

```

### Set working directory
```{r set_wd, include = TRUE}

# SET WORKING DIRECTORY AND LOAD DATA
setwd("C:/Users/samjg/Documents/My_Projects/Pgenerosa_TagSeq_Metabolomics/TagSeq/")
getwd()
path_out = 'C:/Users/samjg/Documents/My_Projects/Pgenerosa_TagSeq_Metabolomics/TagSeq/RAnalysis/DESeq2/output/10cpm'
```



# LOAD DATA: choose ONE of the counts clusters (filtered raw reads or unfiltered raw reads)



### Gene count data - convert to matrix - NOTE ONLY RUN ONE (FILTERED  versus UNFILTERED)
I include both filtered and unfiltered here to compare the results. 
DESEq2 mentions that pre-filtering is not necessary however will reduce the memory/time to run dds and rlog of the models 
Below are both RAW counts filtered and unfiltered to run DESeq2 

# FILTERED COUNTS:  10 CPM (with edgeR) and 50% of samples in each day (days 0, 7, 14, 21)
```{r FILTERED_count.data, include = TRUE}
getwd()
### filtered counts tables - format matrix after upload [from Count_Matrix_Stats.Filter.R] 
d0.counts_matrix  <- read.csv(file="../Data/Filtered_Counts/10cpm_50perc/day0.counts.filtered_10cpm50perc.csv", sep=',', header=TRUE) 
d0.counts_matrix <- data.frame(d0.counts_matrix[,-1], row.names=d0.counts_matrix[,1]) 

d7.counts_matrix  <- read.csv(file="../Data/Filtered_Counts/10cpm_50perc/day7.counts.filtered_10cpm50perc.csv", sep=',', header=TRUE) 
d7.counts_matrix <- data.frame(d7.counts_matrix[,-1], row.names=d7.counts_matrix[,1]) 

d14.counts_matrix <- read.csv(file="../Data/Filtered_Counts/10cpm_50perc/day14.counts.filtered_10cpm50perc.csv", sep=',', header=TRUE) 
d14.counts_matrix <- data.frame(d14.counts_matrix[,-1], row.names=d14.counts_matrix[,1]) 

d21.counts_matrix <- read.csv(file="../Data/Filtered_Counts/10cpm_50perc/day21.counts.filtered_10cpm50perc.csv", sep=',', header=TRUE) 
d21.counts_matrix <- data.frame(d21.counts_matrix[,-1], row.names=d21.counts_matrix[,1]) 


# Pgen annotation master
Pgen_annot <- read.delim(file="../../Seq_details/Panopea-generosa-genes-annotations.txt", header=F) 

```


### load the constant DEGs found up and down in response to primary treatment days 7 14 and 21
This data set is compiled and output in this same markdown file 
cluster == {r Output_tables_and_knitr_summary_tables} (below!)

Objective: use 'Const.DEGs.list' to label callouts in volcano plots 
of the primary treatment effect days 7 14 and 21

```{r load constant DEGs}

# NOTE: this data 
Const.DEGs <- read.csv(file="../Output/DESeq2/10cpm/PrimaryEffect_PersistantDEGs_occurances_10CPM.csv", sep=',', header=TRUE) 
Const.DEGs <- Const.DEGs[,c(2:5)]
names(Const.DEGs)[1] <- "Gene.ID"

annot <- read.delim(file="../../Seq_details/Panopea-generosa-genes-annotations.txt", header=F) 

annot <- annot[,c(1,7:8)]
names(annot) <- c('Gene.ID', 'description', 'GO.terms')

# output the constant DEGs merged with the annotation
merge.constDEGs.annot <- merge(Const.DEGs, annot, by='Gene.ID') # merge with annotation
path_out= 'C:/Users/samjg/Documents/My_Projects/Pgenerosa_TagSeq_Metabolomics/TagSeq/Analysis/Output/DESeq2/10cpm/' # call path
write.csv(merge.constDEGs.annot, paste(path_out,"PrimaryEffect_PersistantDEGs_occurances_10CPM.csv")) # write



# ommit rows containing NA for GO and gene annotation term
Const.DEGs.annotated <- merge.constDEGs.annot[!is.na(merge.constDEGs.annot$GO.terms),]
Const.DEGs.list <- list(Const.DEGs.annotated$Gene.ID) # 12 degs contant upreg or downreg over days 7 14 and 21 WITH Go.term annotation
```


Sample metadata - Experimental treatments/groups
```{r experiment_data, include = TRUE}
### experiment metadata [from Count_Matrix_Stats.Filter.R]  - convert characaters to factors for DESeq2
all.exp_data  <- read.csv(file="../Data/Experiment_Metadata/all.exp.data.csv", sep=',', header=TRUE) %>%   mutate_if(is.character, as.factor)

d0.exp_data   <- read.csv(file="../Data/Experiment_Metadata/day0.exp.data.csv", sep=',', header=TRUE) %>%   mutate_if(is.character, as.factor)

d7.exp_data   <- read.csv(file="../Data/Experiment_Metadata/day7.exp.data.csv", sep=',', header=TRUE) %>%   mutate_if(is.character, as.factor)

d14.exp_data  <- read.csv(file="../Data/Experiment_Metadata/day14.exp.data.csv", sep=',', header=TRUE) %>%   mutate_if(is.character, as.factor)

d21.exp_data  <- read.csv(file="../Data/Experiment_Metadata/day21.exp.data.csv", sep=',', header=TRUE) %>%   mutate_if(is.character, as.factor)
```



# DESeq2 datasets



### DEsign DESeqDatasets
```{r build_dds_datasets}
# ========================================================== 
# DAY 0 FULL MODEL ==  design = ~ Primary_Treatment
# ========================================================== #
#####
d0.metadata <- d0.exp_data %>% dplyr::select(c('Sample.Name', 'Primary_Treatment', 'All_Treatment')) # coondense dataset to build target matrix
d0.metadata <- data.frame(d0.metadata[,-1], row.names=d0.metadata[,1]) # move Sample.Name column as row names  
d0.metadata.mtx <- as.matrix(d0.metadata, row.names="Geoduck.ID") # create matrix 
# check for 'TRUE' in each - check before proceeding  design
d0.metadata.mtx <- d0.metadata.mtx[match(colnames(d0.counts_matrix),rownames(d0.metadata.mtx)), ]
# all(rownames(d0.metadata.mtx) %in% colnames(d0.counts_matrix)) # should be TRUE
# all(rownames(d0.metadata.mtx) == colnames(d0.counts_matrix)) # should be TRUE
# all(rownames(d0.metadata.mtx) == colnames(d0.counts_matrix))  # should be TRUE
# build dds
dds.d0 <- DESeqDataSetFromMatrix(countData = d0.counts_matrix,
                                      colData = d0.metadata.mtx,
                                      design = ~ Primary_Treatment) # DESeq Data Set (dds) - design as ~Primary_Treatment
 
# ========================================================== 
# DAY 7 - FULL MODEL     ==  design = ~ Primary_Treatment+Second_Treament+Primary_Treatment:Second_Treament
#       - ADDITIVE MODEL ==  design = ~ Primary_Treatment+Second_Treament
#       - GROUP    MODEL ==  design = ~ Primary_Treatment+Second_Treament
# ========================================================== #
#####
d7.metadata <- d7.exp_data %>% dplyr::select(c('Sample.Name', 'Primary_Treatment', 'Second_Treament', 'All_Treatment')) # condense dataset to   build target matrix
d7.metadata <- data.frame(d7.metadata[,-1], row.names=d7.metadata[,1]) # move Sample.Name column as row names  
d7.metadata.mtx <- as.matrix(d7.metadata, row.names="Geoduck.ID") # create matrix 
# CKECK THE cts.matrix AND THE exp.data.mtx
d7.metadata.mtx <- d7.metadata.mtx[match(colnames(d7.counts_matrix),rownames(d7.metadata.mtx)), ]
# all(rownames(d7.metadata.mtx) %in% colnames(d7.counts_matrix)) # should be TRUE
# all(rownames(d7.metadata.mtx) == colnames(d7.counts_matrix)) # should be TRUE
# all(rownames(d7.metadata.mtx) == colnames(d7.counts_matrix))  # should be TRUE
# build dds
# FULL MODEL 
dds.d7 <- DESeqDataSetFromMatrix(countData = d7.counts_matrix,
                                  colData = d7.metadata.mtx,
                                  design = ~ Second_Treament+Primary_Treatment+Second_Treament:Primary_Treatment) # DESeq Data Set (dds) - design as ~Primary_Treatment

# ADDITIVE MODEL  
dds.d7.main <- DESeqDataSetFromMatrix(countData = d7.counts_matrix,
                                       colData = d7.metadata.mtx,
                                       design = ~ Second_Treament+Primary_Treatment)
# GROUP MODEL- note! ~group -1` removed the intercept to call each pariwise interaction
dds.d7.group <- DESeqDataSetFromMatrix(countData = d7.counts_matrix,
                                        colData = d7.metadata.mtx,
                                        design = ~ All_Treatment-1) 
# ========================================================== 
# DAY 14  - FULL MODEL     ==  design = ~ Primary_Treatment+Second_Treament+Primary_Treatment:Second_Treament
#       - ADDITIVE MODEL ==  design = ~ Primary_Treatment+Second_Treament
#       - GROUP    MODEL ==  design = ~ Primary_Treatment+Second_Treament
# ========================================================== #
#####
d14.metadata <- d14.exp_data %>% dplyr::select(c('Sample.Name', 'Primary_Treatment', 'Second_Treament', 'All_Treatment')) # condense dataset to build target matrix
d14.metadata <- data.frame(d14.metadata[,-1], row.names=d14.metadata[,1]) # move Sample.Name column as row names  
d14.metadata.mtx <- as.matrix(d14.metadata, row.names="Geoduck.ID") # create matrix 
# CKECK THE cts.matrix AND THE exp.data.mtx
d14.metadata.mtx <- d14.metadata.mtx[match(colnames(d14.counts_matrix),rownames(d14.metadata.mtx)), ]
# all(rownames(d14.metadata.mtx) %in% colnames(d14.counts_matrix)) # should be TRUE
# all(rownames(d14.metadata.mtx) == colnames(d14.counts_matrix)) # should be TRUE
# all(rownames(d14.metadata.mtx) == colnames(d14.counts_matrix))  # should be TRUE
# build dds
# FULL MODEL 
dds.d14 <- DESeqDataSetFromMatrix(countData = d14.counts_matrix,
                                     colData = d14.metadata.mtx,
                                     design = ~ Second_Treament+Primary_Treatment+Second_Treament:Primary_Treatment) # DESeq Data Set (dds) - design as ~Primary_Treatment
# GROUP MDOEL
dds.d14.group <- DESeqDataSetFromMatrix(countData = d14.counts_matrix,
                                        colData = d14.metadata.mtx,
                                        design = ~ All_Treatment-1) 
# ADDITIVE MODEL 
dds.d14.main <- DESeqDataSetFromMatrix(countData = d14.counts_matrix,
                                       colData = d14.metadata.mtx,
                                       design = ~ Second_Treament+Primary_Treatment) 
# ========================================================== 
# DAY 21 -  FULL MODEL == design = ~Primary_Treatment+Second_Treament+Third_Treatment+Primary_Treatment:Second_Treament+Primary_Treatment:Third_Treatment+Second_Treament:Third_Treatment
#---- INPUT: Experiment/design matrix == 'exp.data.d21'
#---- INPUT: TagSeq count matrix == 'cts.matrix.d21.filtered' (3 CPM in 50% samples)
# ========================================================== 
#####
d21.metadata <- d21.exp_data %>% dplyr::select(c('Sample.Name', 'Primary_Treatment','Second_Treament', 'Third_Treatment', 'All_Treatment')) # coondense dataset to build target matrix
d21.metadata <- data.frame(d21.metadata[,-1], row.names=d21.metadata[,1]) # move Sample.Name column as row names  
d21.metadata.mtx <- as.matrix(d21.metadata, row.names="Geoduck.ID") # create matrix 
# CKECK THE cts.matrix AND THE exp.data.mtx
d21.metadata.mtx <- d21.metadata.mtx[match(colnames(d21.counts_matrix),rownames(d21.metadata.mtx)), ]
# all(rownames(d21.metadata.mtx) %in% colnames(d21.counts_matrix)) # should be TRUE
# all(rownames(d21.metadata.mtx) == colnames(d21.counts_matrix)) # should be TRUE
# all(rownames(d21.metadata.mtx) == colnames(d21.counts_matrix))  # should be TRUE
# build dds
# FULL MODEL
dds.d21 <- DESeqDataSetFromMatrix(countData = d21.counts_matrix,
                                        colData = d21.metadata.mtx,
                                        design = ~ Primary_Treatment+
                                        Second_Treament+
                                        Third_Treatment+
                                        Primary_Treatment:Second_Treament+
                                        Primary_Treatment:Third_Treatment+
                                        Second_Treament:Third_Treatment) 
# GROUP MDOEL
dds.d21.group <- DESeqDataSetFromMatrix(countData = d21.counts_matrix,
                                        colData = d21.metadata.mtx,
                                        design = ~ All_Treatment-1) 
# ADDITIVE MODEL 
dds.d21.main <- DESeqDataSetFromMatrix(countData = d21.counts_matrix,
                                       colData = d21.metadata.mtx,
                                       design = ~ Third_Treatment+Second_Treament+Primary_Treatment) 

```

### Plot reads per gene and reads per sample in each dds
```{r ddsPlots_read.gene_reads.sample}
# DAY 0 PLOT THE (1) READS PER SAMPLE (2) READS PER GENE - FOR EACH DDS
d0.nsamples <- ncol(counts(dds.d0)) # Number of samples - for the plot label
d0.rps <- qplot(colSums(counts(dds.d0))) +# reads of reads per sample 
  labs(x = "Mapped reads per sample", y = "Number of samples",
       title = "Day 0: Mapped reads per sample") +
  geom_label(aes(x = 8e5, y = 1, label = paste(d0.nsamples, "samples")))
d0.ngenes <- nrow(counts(dds.d0)) # Number of genes
d0.ngenes_min <-min(rowSums(counts(dds.d0))) #  minimum reads
d0.ngenes_mean.min <- min(rowMeans2(counts(dds.d0))) # minimum row mean
d0.ngenes_max <-max(rowSums(counts(dds.d0))) #  maximum reads
d0.rpg <- qplot(log10(rowSums(counts(dds.d0))), bins = 16) + # number of reads per gene
  labs(x = "log10(Mapped reads per gene)", y = "Number of genes",
       title = "Day 0: Mapped reads per gene") +
  geom_label(aes(x = 4, y = 3000, label = paste(d0.ngenes, "total genes"))) +
  geom_label(aes(x = 4, y = 2700, label = paste("max count =", d0.ngenes_max))) +
  geom_label(aes(x = 4, y = 2400, label = paste("min count =", d0.ngenes_min))) +
  geom_label(aes(x = 4, y = 2100, label = paste("min mean =", d0.ngenes_mean.min)))
d0.countfig <- plot_grid(d0.rps, d0.rpg)

# DAY 7 PLOT THE (1) READS PER SAMPLE (2) READS PER GENE - FOR EACH DDS
d7.nsamples <- ncol(counts(dds.d7)) # Number of samples - for the plot label
d7.rps <- qplot(colSums(counts(dds.d7))) +# reads of reads per sample 
  labs(x = "Mapped reads per sample", y = "Number of samples",
       title = "Day 7: Mapped reads per sample") +
  geom_label(aes(x = 6e5, y = 4, label = paste(d7.nsamples, "samples")))
d7.ngenes <- nrow(counts(dds.d7)) # Number of genes
d7.ngenes_min <-min(rowSums(counts(dds.d7))) #  minimum reads
d7.ngenes_max <-max(rowSums(counts(dds.d7))) #  maximum reads
d7.ngenes_mean.min <- min(rowMeans2(counts(dds.d7))) # minimum row mean
d7.rpg <- qplot(log10(rowSums(counts(dds.d7))), bins = 16) + # number of reads per gene
  labs(x = "log10(Mapped reads per gene)", y = "Number of genes",
       title = "Day 7: Mapped reads per gene") +
  geom_label(aes(x = 4, y = 3000, label = paste(d7.ngenes, "total genes"))) +
  geom_label(aes(x = 4, y = 2700, label = paste("max count =", d7.ngenes_max))) +
  geom_label(aes(x = 4, y = 2400, label = paste("min count =", d7.ngenes_min)))  +
  geom_label(aes(x = 4, y = 2100, label = paste("min mean =", d7.ngenes_mean.min)))
d7.countfig <- plot_grid(d7.rps, d7.rpg)


# DAY 14 PLOT THE (1) READS PER SAMPLE (2) READS PER GENE - FOR EACH DDS
d14.nsamples <- ncol(counts(dds.d14)) # Number of samples - for the plot label
d14.rps <- qplot(colSums(counts(dds.d14))) +# reads of reads per sample 
  labs(x = "Mapped reads per sample", y = "Number of samples",
       title = "Day 14: Mapped reads per sample") +
  geom_label(aes(x = 6e5, y = 4, label = paste(d14.nsamples, "samples")))
d14.ngenes <- nrow(counts(dds.d14)) # Number of genes
d14.ngenes_min <-min(rowSums(counts(dds.d14))) #  minimum reads
d14.ngenes_max <-max(rowSums(counts(dds.d14))) #  maximum reads
d14.ngenes_mean.min <- min(rowMeans2(counts(dds.d14))) # minimum row mean
d14.rpg <- qplot(log10(rowSums(counts(dds.d14))), bins = 16) + # number of reads per gene
  labs(x = "log10(Mapped reads per gene)", y = "Number of genes",
       title = "Day 14: Mapped reads per gene") +
  geom_label(aes(x = 4, y = 3000, label = paste(d14.ngenes, "total genes"))) +
  geom_label(aes(x = 4, y = 2700, label = paste("max count =", d14.ngenes_max))) +
  geom_label(aes(x = 4, y = 2400, label = paste("min count =", d14.ngenes_min))) +
  geom_label(aes(x = 4, y = 2100, label = paste("min mean =", d14.ngenes_mean.min)))
d14.countfig <- plot_grid(d14.rps, d14.rpg)

# DAY 21 PLOT THE (1) READS PER SAMPLE (2) READS PER GENE - FOR EACH DDS
d21.nsamples <- ncol(counts(dds.d21)) # Number of samples - for the plot label
d21.rps <- qplot(colSums(counts(dds.d21))) +# reads of reads per sample 
  labs(x = "Mapped reads per sample", y = "Number of samples",
       title = "Day 21: Mapped reads per sample") +
  geom_label(aes(x = 6e5, y = 11, label = paste(d21.nsamples, "samples")))
d21.ngenes <- nrow(counts(dds.d21)) # Number of genes
d21.ngenes_min <-min(rowSums(counts(dds.d21))) #  minimum reads
d21.ngenes_max <-max(rowSums(counts(dds.d21))) #  maximum reads
d21.ngenes_mean.min <- min(rowMeans2(counts(dds.d21))) # minimum row mean
d21.countfig <- qplot(log10(rowSums(counts(dds.d21))), bins = 16) + # number of reads per gene
  labs(x = "log10(Mapped reads per gene)", y = "Number of genes",
       title = "Day 21: Mapped reads per gene") +
  geom_label(aes(x = 4, y = 3000, label = paste(d21.ngenes, "total genes"))) +
  geom_label(aes(x = 4, y = 2700, label = paste("max count =", d21.ngenes_max))) +
  geom_label(aes(x = 4, y = 2400, label = paste("min mean =", d21.ngenes_min))) +
  geom_label(aes(x = 4, y = 2100, label = paste("min mean =", d21.ngenes_mean.min)))
d21.countfig <- plot_grid(d21.rps, d21.countfig)

png("../Data/DESeq2/dds_GeneCounts_10cpm.png", 1000, 1000, pointsize=20)
plot_grid(d0.countfig, d7.countfig, d14.countfig, d21.countfig)
dev.off()
```

### RUN DESEQ2 model - wait for this to complete...
# NOTE: this will be longer dependent on filtered vs. unfiltered raw read count data for dds objects 
```{r run_dds}

dds.d0 <- DESeq(dds.d0) # wait for this to complete....

# RUN DESEQ2 model - view all the pariwise comparisons
dds.d7 <- DESeq(dds.d7) #  full model                 wait for this to complete....
dds.d7.group <- DESeq(dds.d7.group) # group model     wait for this to complete....
dds.d7.main <- DESeq(dds.d7.main) # main effect model    

# RUN DESEQ2 model - view all the pariwise comparisons
dds.d14 <- DESeq(dds.d14) #  full model               wait for this to complete....
dds.d14.group <- DESeq(dds.d14.group) # group model   wait for this to complete....
dds.d14.main <- DESeq(dds.d14.main) # main effect model    

# RUN DESEQ2 model - view all the pariwise comparisons
dds.d21 <- DESeq(dds.d21) #  full model                  wait for this to complete....
dds.d21.group <- DESeq(dds.d21.group) # group model      wait for this to complete....
dds.d21.main <- DESeq(dds.d21.main) # main effect model    

```

# Analysis (DEG tables)
Rationale for thresholds in the following cluster:
With the thousands of pairwise comparisons between genes, it is critical to acknowledge the prevalence of spurious discoveries that must be guarded against!
In the following, I will run each model with 'alpha = 0.05' to assume a False Discovery Rate of 5%. I reinforce this thinking by graphing all pvalues 
as a histogram to show the distrubution - abline(s) show the FDR cut-off and the true padj values at 0.05 and 0.1 thresholds. 

Strader et al. 2018 - 48 total samples, genes with mean count <10 across all samples werre ommitted
                    - rlog transformated (rlog fxn)
                    - adonic and vegan package to assess multivariate analysis of variance 
                    - Used a 10% false discovery rate threshold
                    - GO enrichments performend using the stat value output from DESeq2 using Mann Whitney U test (here: https://github.com/z0on/ GO_MWU)
Johnson and Kelly 2021 - 40 samples 
                       - used edgeR to assess pairwise changes in gene expression between treatments
                       - filtered genes with fewer than 3 CPM across 50%  *n=20_ of all samples
                       used R program vegan to assess PCoA with Euclidean distances from log2+1 transformed normalized counts from cpm function in egde R
                       - used Mann Whitney U test to identify enriched ontologies - than this also after identifying modules via WGCNA with physiological measurements


NOTE: 
Read 'Interactions' here 
https://bioconductor.riken.jp/packages/3.8/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#interactions
When adding the interaction term x:y to a dds model, the main condition effect only represents the effect of that ccondition
for the REFERENCe level in the other. Example: design = ~primary + second + primary:second
levels = A M (primary) & A  M S (second)
The main effect of primary as Primary_A_v_M in this model is ONLY in response to A (ref level) for the treatment 'second' 




# Main Treatment effects: padj < 0.05; FDR 5%; LFC >1 (<-1)
Primary: Ambient v Moderate; 
Second:  Ambient v Moderate,  Ambient v Severe, & Moderate v Severe; 
Third:   Ambient v Moderate

#Day 0 - NOTE: this model is run is the oppositoe direction of the other (M v A as opposed to A v M)
```{r Day0_Main_effects}

# DAY 0 - main model
# resd0.primary <- results(dds.d0, name="Primary_Treatment_M_vs_A", alpha = 0.05)
resd0.primary <- results(dds.d0, contrast=c("Primary_Treatment", "M", "A"), alpha = 0.05)
d0.numDEGs_padj <-data.frame(table(resd0.primary$padj<0.05))[2,2] # DEGs - NOT considering LFC - just p adj
resd0.primary.ordered <- resd0.primary[order(resd0.primary$padj), ] # Order by adjusted p-value
d0.num.DownReg <- sum((resd0.primary.ordered$log2FoldChange[1:d0.numDEGs_padj] > 0) == TRUE) #  LFC >= 1
d0.num.UpReg <- sum((resd0.primary.ordered$log2FoldChange[1:d0.numDEGs_padj] < 0) == TRUE) # LFC >= 1
d0.total <- sum(d0.num.DownReg,d0.num.UpReg) # sum of DEGs with the criteria pdj < 0.05 + LFC>1 (< -1)
# histogram of FDR
png("../Output/DESeq2/10cpm/Day0/Day0.FDR_hist.png", 1000, 1000, pointsize=20)
hist(resd0.primary$pvalue, breaks=20, col="grey") # view histogram
abline(h=c( (nrow(resd0.primary)*0.05), 
            ((table(resd0.primary$padj < 0.1)[2]) + (nrow(resd0.primary)*0.1)),
            ((table(resd0.primary$padj < 0.05)[2]) + (nrow(resd0.primary)*0.05)) ),
                  col=c("blue", "red", "red"), lty=c(1,2,1), lwd=c(1, 3, 1)) # add line at 
dev.off()


# Write results - covert to as.data.frame for the ordered results
resdata.d0.primary  <- merge(as.data.frame(resd0.primary.ordered), as.data.frame(counts(dds.d0, normalized=TRUE)), by="row.names", sort=FALSE) ## Merge with normalized count data
names(resdata.d0.primary)[1] <- "Gene"
resdata.d0.primary <- resdata.d0.primary[order(resdata.d0.primary$padj), ] # Order by adjusted p-value
path_out.d0 = 'C:/Users/samjg/Documents/My_Projects/Pgenerosa_TagSeq_Metabolomics/TagSeq/Analysis/Output/DESeq2/10cpm/Day0/' # call path
write.csv(resdata.d0.primary, paste(path_out.d0,"Day0.Primary_DESeq2results.csv")) # write


#MASTER MAIN EFFECTS TABLE FOR OUTPUT  ======================================================================================================== #
# rbind all significant DESeq results in main effect models 
Day0_MainEffectMaster <- (as.data.frame(resdata.d0.primary)[c(1,3,7)]   %>% dplyr::mutate(Day_treatment = "Day0", Treatment = "Primary", Pairwise = "AvM") %>%  dplyr::filter(padj < 0.05) )
# load annotation to merge by GeneID 
Pgen_annot <- read.delim(file="../../Seq_details/Panopea-generosa-genes-annotations.txt", header=F) 
Pgen_annot <- Pgen_annot[,c(1,5,7)]
names(Pgen_annot) <- c('Gene', 'Uniprot', 'Gene_terms_EC')
# merge the master with the Pgen_annot and save the csv file 
Day0_MainEffectMaster_2 <- merge(Day0_MainEffectMaster, Pgen_annot, by ="Gene")
path_out.d0 = 'C:/Users/samjg/Documents/My_Projects/Pgenerosa_TagSeq_Metabolomics/TagSeq/Analysis/Output/DESeq2/10cpm/Day0/' # call path to save
write.csv(Day0_MainEffectMaster_2, paste(path_out.d0,"Day0.AllPrimaryDEGs_DESeq2results.csv")) # write



# volcano plot 
png("../Output/DESeq2/10cpm/Day0/Day0.PrimaryTreatment-VolcanoPlot.png", 1000, 1000, pointsize=20)
EnhancedVolcano(resd0.primary,
                lab = rownames(resd0.primary),
                x = 'log2FoldChange',
                y = 'padj',
                title = 'Day 0: Primary Treatment (M v A)',
                subtitle = "DESeq2 - Differential expression",
                FCcutoff = 1,
                pCutoff = 0.05,
                pointSize = 4.0,
                labSize = 6.0,
                colAlpha = 1,
                legendPosition = 'right',
                legendLabSize = 12,
                legendIconSize = 4.0,
                widthConnectors = 0.75)
dev.off()


# Plot dispersions 
png("../Output/DESeq2/10cpm/Day0/Day0.dispersions.png", 1000, 1000, pointsize=20)
plotDispEsts(dds.d0, main="Day 0 dispersions")
dev.off()

# Data transformations for heatmap and PCA visuals ======================================================================================= #
# rlog - regularized log transformation of origin count data to log2 scale - fit for each sample and dist. of coefficients in the data
rlog.d0<- rlogTransformation(dds.d0) # rlog transform (regularized log)

# diagnostiscs
png("../Output/DESeq2/10cpm/Day0/Day0.rlog_histogram.png", 1000, 1000, pointsize=20)# diagnostics of transformation # Histogram and sd plot
hist(assay(rlog.d0)) # view histogram 
dev.off()
png("../Output/DESeq2/10cpm/Day0/Day0.rlog_mean_sd.png", 1000, 1000, pointsize=20)
meanSdPlot(assay(rlog.d0)) # shows the sd y axis (sq root of varaince in all samples) - flat curve may seem like a goals, BUT may be unreasonable in cases with MANY true DEGs from experimental conditions
dev.off()

# PCA plot rlog ------------------------ #
pcaData_d0 <- plotPCA(rlog.d0, intgroup = "Primary_Treatment", returnData = TRUE)
percentVar <- round(100 * attr(pcaData_d0, "percentVar"))
png("../Output/DESeq2/10cpm/Day0/Day0.rlog_PCA.png", 1000, 1000, pointsize=20)
ggplot(pcaData_d0, aes(x = PC1, y = PC2, color = Primary_Treatment, label=name)) +
  scale_shape_manual(values = c(4, 19, 17)) +
  geom_text(aes(label=name),hjust=0.2, vjust=1.4, size=5) +
  geom_point(size =6) +
  theme_classic() +
  theme(text = element_text(size=15)) +
  ggtitle("PCA: Day0 (rlog)") +
  xlab(paste0("PC1: ", percentVar[1], "% variance")) +
  ylab(paste0("PC2: ", percentVar[2], "% variance")) +
  coord_fixed()
dev.off()


# Plot pheatmap map rlog------------------------ #
# save pheatmap fxn
save_pheatmap <- function(x, filename, width=1000, height=960) { # template for saving pheatmap outputs
  stopifnot(!missing(x))
  stopifnot(!missing(filename))
  png(filename,width = width, height=height)
  grid::grid.newpage()
  grid::grid.draw(x$gtable)
  dev.off()
}
# parameters for the pheatmap
topgenes.IDs <- head(rownames(resd0.primary.ordered),14) # call the first 13 genes that are differentially expressed, FDR < 0.05 correct pdj
topgenes_rlog.counts <- assay(rlog.d0)[topgenes.IDs,] 
topgenes_corrected <- topgenes_rlog.counts - rowMeans(topgenes_rlog.counts) # substract from the row mean to get +/- 0 to normalizze and ease the asthetic
df_annot.col <- as.data.frame(colData(dds.d0)[c("Primary_Treatment")])


# colorblind palette
# colorBlindBlack8  <- c("#000000", "#E69F00", "#56B4E9", "#009E73", 
#                        "#F0E442", "#0072B2", "#D55E00", "#CC79A7")
# pie(rep(1, 8), col = colorBlindBlack8)

ann_colors = list(
    Primary_Treatment = c(A="turquoise", M="coral1"))

d0.rlog.heatmap <- pheatmap(topgenes_corrected, 
                            annotation_col=df_annot.col, 
                            main = "Day 0  Ambient versus Moderate (14 total DEGs)",
                            cutree_cols = 2,
                            cutree_rows = 2,
                            annotation_legend = TRUE,
                            annotation_colors = ann_colors,
                            show_rownames = FALSE,
                            labels_col=df_annot.col$Primary_Treatment, 
                            angle_col = 0,
                            fontsize = 8,
                            legend = TRUE)
#Save ppheatmap
save_pheatmap(d0.rlog.heatmap, filename = "../Output/DESeq2/10cpm/Day0/Day0.PRIMARY.rlog_heatmap.png")

```

#Day 7
```{r Day7_Main_effects}
# TEST! (read below)
# 'Interaction' term tests - run this to see that the 'Primary Treatment' main effect in NOT represented in the full effect model
test1 <- results(dds.d7, contrast = c("Primary_Treatment", "M", "A"),  alpha= 0.05) # D/N addresses Primary M vs A in reference level Secondary!!! 
# thus the model above measures the pairwsie comparisons of MA vs. AA
table(test1$padj<0.05) #  10 DEGs
test2 <- results(dds.d7.group, contrast = list(("All_TreatmentMA"),("All_TreatmentAA")),  alpha= 0.05) # look here at the group - same genes are DEGs
table(test2$padj<0.05) #  10 DEGs (same as above)


# Primary Treatmetn A v M ========================================================================================================== #



# group by contrast
res.d7.primary <- results(dds.d7.group, contrast = list(c("All_TreatmentAA", "All_TreatmentAM", "All_TreatmentAS"),c("All_TreatmentMA", "All_TreatmentMM", "All_TreatmentMS")),  alpha= 0.05)
table(res.d7.primary$padj<0.05) #  108 DEGs 

# main effect model (mo interaction term)
res.d7.primary_main <- results(dds.d7.main, contrast = c("Primary_Treatment", "A", "M"),  alpha= 0.05)
table(res.d7.primary_main$padj<0.05) #  98 DEGs 

# full model 
res.d7.primary_fullmod <- results(dds.d7, contrast = c("Primary_Treatment", "A", "M"),  alpha= 0.05)
table(res.d7.primary_fullmod$padj<0.05) #  10 DEGs 

# Group contrast - this tests the full model effect above - appears the interaction term in fullmod is the pariwise effect under the reference level in treatment 2
res.d7.primary_grouptest<- results(dds.d7.group, contrast = list(("All_TreatmentAA"),("All_TreatmentMA")),  alpha= 0.05) # look here at the group - same genes are DEGs
table(res.d7.primary_grouptest$padj<0.05) #  10 DEGs (same as above?)



# Secondary treatment A v M ========================================================================================================= #



# group by contrast
res.d7.Second_AM <- results(dds.d7.group, contrast = list(c("All_TreatmentAA", "All_TreatmentMA"),c("All_TreatmentAM", "All_TreatmentMM")),  alpha= 0.05)
table(res.d7.Second_AM$padj<0.05) # 106 DEGs

# main effect model (mo interaction term)
res.d7.Second_AM_main <- results(dds.d7.main, contrast = c("Second_Treament", "A", "M"),  alpha= 0.05)
table(res.d7.Second_AM_main$padj<0.05) #  94 DEGs 

# full model 
res.d7.Second_AM_fullmod <- results(dds.d7, contrast = c("Second_Treament", "A", "M"),  alpha= 0.05) # appears this is actually the pairwise comparisons between AA and AM
table(res.d7.Second_AM_fullmod$padj<0.05) #  0 DEGs 

# Group contrast - test the putative pariwise comparison with the ref level in fullmod
res.d7.Second_AM_grouptest<- results(dds.d7.group, contrast = list(("All_TreatmentAA"),("All_TreatmentAM")),  alpha= 0.05) # look here at the group constrasts
table(res.d7.Second_AM_grouptest$padj<0.05) #  0 DEGs (same as above?)


# after running all the group interaction slater in this script... I found that the majority of the A v M second treatment effect 
# appears to be ddriven by AM  vs. MA treatment groups. Here is the table supporting this.
res.d7.Second_MA_AM <- results(dds.d7.group, contrast = list(c("All_TreatmentMA"),c("All_TreatmentAM")),  alpha= 0.05)
table(res.d7.Second_MA_AM$padj<0.05) # 168 DEGs - 168 total DEGs in this pairwise group model! 
# what does this say?? The second treatment had an effect on differneital expression, and the weight of this effect 
# is driven by the recovery of the moderate primed animals versus the initial encounter of the ambient (naive) animals! 
# outpt this table to Day 7 folder to run a venn diagram agains the second ambient vs. moderate (106 DEGs total) to look for overlap here


# Secondary treatment A v S ========================================================================================================= #



res.d7.Second_AS <- results(dds.d7.group, contrast = list(c("All_TreatmentAA", "All_TreatmentMA"),c("All_TreatmentAS", "All_TreatmentMS")),  alpha= 0.05)
table(res.d7.Second_AS$padj<0.05) #  2 DEGs

# main effect model (mo interaction term)
res.d7.Second_AS_main <- results(dds.d7.main, contrast = c("Second_Treament", "A", "S"),  alpha= 0.05)
table(res.d7.Second_AS_main$padj<0.05) #  2 DEGs 

# full model 
res.d7.Second_AS_fullmod <- results(dds.d7, contrast = c("Second_Treament", "A", "S"),  alpha= 0.05) # appears this is actually the pairwise comparisons between AA and AS
table(res.d7.Second_AS_fullmod$padj<0.05) #  10 DEGs 

# Group contrast - test the putative pariwise comparison with the ref level in fullmod
res.d7.Second_AS_grouptest<- results(dds.d7.group, contrast = list(("All_TreatmentAA"),("All_TreatmentAS")),  alpha= 0.05) # look here at the group constrasts
table(res.d7.Second_AS_grouptest$padj<0.05) #  10 DEGs (same as above?)



# Secondary treatment M v S ========================================================================================================= #



res.d7.Second_MS <- results(dds.d7.group, contrast = list(c("All_TreatmentAM", "All_TreatmentMM"),c("All_TreatmentAS", "All_TreatmentMS")),  alpha= 0.05)
table(res.d7.Second_MS$padj<0.05) #  13 DEGs

# main effect model (mo interaction term)
res.d7.Second_MS_main <- results(dds.d7.main, contrast = c("Second_Treament", "M", "S"),  alpha= 0.05)
table(res.d7.Second_MS_main$padj<0.05) #  14 DEGs 

# full model 
res.d7.Second_MS_fullmod <- results(dds.d7, contrast = c("Second_Treament", "M", "S"),  alpha= 0.05) # appears this is actually the pairwise comparisons between AM and AS
table(res.d7.Second_MS_fullmod$padj<0.05) #  0 DEGs 

# Group contrast - test the putative pariwise comparison with the ref level in fullmod
res.d7.Second_MS_grouptest<- results(dds.d7.group, contrast = list(("All_TreatmentAM"),("All_TreatmentAS")),  alpha= 0.05) # look here at the group constrasts
table(res.d7.Second_MS_grouptest$padj<0.05) #  0 DEGs (same as above?)





# TABLES  &  Volcano Plots








# Primary Treatmetn M v A ========================================================================================================= #

# DEGs total for the primary treatment effect
d7.numDEGs_padj               <- data.frame(table(res.d7.primary$padj<0.05))[2,2] # DEGs - NOT considering LFC - just p adj
res.d7_primary.ordered <- res.d7.primary[order(res.d7.primary$padj), ] # Order by adjusted p-value - 110 DEGs
d7.num.UpReg                  <- sum((res.d7_primary.ordered$log2FoldChange[1:d7.numDEGs_padj] > 0) == TRUE) #  LFC >= 1
d7.num.DownReg                <- sum((res.d7_primary.ordered$log2FoldChange[1:d7.numDEGs_padj] < 0) == TRUE) # LFC >= 1
d7.total                      <- sum(d7.num.DownReg,d7.num.UpReg) # 108 sum of DEGs with the criteria pdj < 0.05 + LFC>1 (< -1)

# Histogram(s) 
png("../Output/DESeq2/10cpm/Day7/Day7.Primary.FDR_hist.png", 1000, 1000, pointsize=20)
hist(res.d7.primary$pvalue, breaks=20, col="grey") # view histogram,  
abline(h=c( (nrow(res.d7.primary)*0.05), 
            ((table(res.d7.primary$padj < 0.1)[2]) + (nrow(res.d7.primary)*0.05)),
            ((table(res.d7.primary$padj < 0.05)[2]) + (nrow(res.d7.primary)*0.05)) ),
                  col=c("blue", "red", "red"), lty=c(1,2), lwd=c(1, 3)) # add line at 
dev.off()

# Write results csv ( here I output the four tests of the primary treatment to run a Venn diagram on shared genes - sanity check to test what these models are doing differently)
# group results 
resdata.d7.primary  <- merge(as.data.frame(res.d7_primary.ordered), as.data.frame(counts(dds.d7, normalized=TRUE)), by="row.names", sort=FALSE) ## Merge with normalized count data
names(resdata.d7.primary)[1] <- "Gene"
resdata.d7.primary <- resdata.d7.primary[order(resdata.d7.primary$padj), ] # Order by adjusted p-value
path_out.d7 = 'C:/Users/samjg/Documents/My_Projects/Pgenerosa_TagSeq_Metabolomics/TagSeq/Analysis/Output/DESeq2/10cpm/Day7/' # call path
write.csv(resdata.d7.primary, paste(path_out.d7,"Day7.PrimaryDEGs_AvM__DESeq2results.csv")) # write

# main effect model
resdata.d7.primary_main  <- merge(as.data.frame(res.d7.primary_main), as.data.frame(counts(dds.d7, normalized=TRUE)), by="row.names", sort=FALSE) ## Merge with normalized count data
names(resdata.d7.primary_main)[1] <- "Gene"
resdata.d7.primary_main <- resdata.d7.primary_main[order(resdata.d7.primary_main$padj), ] # Order by adjusted p-value
path_out.d7.tests = 'C:/Users/samjg/Documents/My_Projects/Pgenerosa_TagSeq_Metabolomics/TagSeq/Analysis/Output/DESeq2/10cpm/Day7/D7_res_model_tests/' # call path
write.csv(resdata.d7.primary_main, paste(path_out.d7.tests,"Day7.Primary_res_main.mod.csv")) # write

# full model
resdata.d7.primary_fullmod  <- merge(as.data.frame(res.d7.primary_fullmod), as.data.frame(counts(dds.d7, normalized=TRUE)), by="row.names", sort=FALSE) ## Merge with normalized count data
names(resdata.d7.primary_fullmod)[1] <- "Gene"
resdata.d7.primary_fullmod <- resdata.d7.primary_fullmod[order(resdata.d7.primary_fullmod$padj), ] # Order by adjusted p-value
path_out.d7.tests = 'C:/Users/samjg/Documents/My_Projects/Pgenerosa_TagSeq_Metabolomics/TagSeq/Analysis/Output/DESeq2/10cpm/Day7/D7_res_model_tests/' # call path
write.csv(resdata.d7.primary_fullmod, paste(path_out.d7.tests,"Day7.Primary_res_full.mod.csv")) # write

# group test of full model
resdata.d7.primary_grouptest  <- merge(as.data.frame(res.d7.primary_grouptest), as.data.frame(counts(dds.d7, normalized=TRUE)), by="row.names", sort=FALSE) ## Merge with normalized count data
names(resdata.d7.primary_grouptest)[1] <- "Gene"
resdata.d7.primary_grouptest <- resdata.d7.primary_grouptest[order(resdata.d7.primary_grouptest$padj), ] # Order by adjusted p-value
path_out.d7.tests = 'C:/Users/samjg/Documents/My_Projects/Pgenerosa_TagSeq_Metabolomics/TagSeq/Analysis/Output/DESeq2/10cpm/Day7/D7_res_model_tests/' # call path
write.csv(resdata.d7.primary_grouptest, paste(path_out.d7.tests,"Day7.Primary_res_grouptest.mod.csv")) # write

# volcano plot 
# Set max.overlapos to infinite globally to allow muliple labels that overlap in the volcano plots 
options(ggrepel.max.overlaps = Inf) # from library(ggrepel)
# Const.DEGs.list # copy this list for selectLab callouts
png("../Output/DESeq2/10cpm/Day7/Day7.PrimaryTreatment-VolcanoPlot.png", 1000, 1000, pointsize=20)
EnhancedVolcano(res.d7.primary,
                            lab = row.names(res.d7.primary),
                            x = 'log2FoldChange',
                            y = 'padj',
                            title = 'Day 7: Primary Treatment (A v M)',
                            subtitle = "DESeq2 - Differential expression",
                            selectLab = c("PGEN_.00g057310", "PGEN_.00g279590", "PGEN_.00g092270", "PGEN_.00g092270", "PGEN_.00g266130"), # call persistant DEGs (day 7,14,21) with mean LFC > 5
                            # selectLab = c("PGEN_.00g003940", "PGEN_.00g057310", "PGEN_.00g092270" ,
                            #               "PGEN_.00g169480", "PGEN_.00g188250", "PGEN_.00g199810",
                            #               "PGEN_.00g211250" ,"PGEN_.00g266130" ,"PGEN_.00g272710",
                            #               "PGEN_.00g279590", "PGEN_.00g348220",
                            #               "PGEN_.00g348570"), # the list of genes from 'Const.DEGs.list'
                            xlab = bquote(~Log[2]~ 'fold change'),
                            ylab = bquote(~-Log[10]~ 'Padj'),
                            pCutoff = 0.05,
                            FCcutoff = 1.0,
                            pointSize = 3.0,
                            labSize = 3.0,
                            labCol = 'black',
                            labFace = 'bold',
                            boxedLabels = TRUE,
                          # parseLabels = TRUE,
                            col = c('black', '#E69F00', '#D55E00', '#0072B2'), # colorblindness pallette
                            colAlpha = 1,
                            legendPosition = 'right',
                            legendLabSize = 10,
                            legendIconSize = 3.0,
                            drawConnectors = TRUE,
                            widthConnectors = 0.2,
                            colConnectors = 'black')
dev.off()







# Secondary treatment A v M ========================================================================================================= #

# DEGs total for the primary treatment effect
d7.numDEGs_padj_S.AM               <- data.frame(table(res.d7.Second_AM$padj<0.05))[2,2] # DEGs - NOT considering LFC - just p adj
res.d7.Second_AM_ordered           <- res.d7.Second_AM[order(res.d7.Second_AM$padj), ] # Order by adjusted p-value - 110 DEGs
d7.num.UpReg_S.AM                  <- sum((res.d7.Second_AM_ordered$log2FoldChange[1:d7.numDEGs_padj_S.AM] > 0) == TRUE) #  LFC >= 1
d7.num.DownReg_S.AM                <- sum((res.d7.Second_AM_ordered$log2FoldChange[1:d7.numDEGs_padj_S.AM] < 0) == TRUE) # LFC >= 1
d7.total_Sec.AM                    <- sum(d7.num.UpReg_S.AM,d7.num.DownReg_S.AM) # 106 TOTAL DEGs with the criteria pdj < 0.05 + LFC>1 (< -1)

# Histogram(s)
png("../Output/DESeq2/10cpm/Day7/Day7.Second_AM.FDR_hist.png", 1000, 1000, pointsize=20)
hist(res.d7.Second_AM$pvalue, breaks=20, col="grey") # view histogram,  
abline(h=c( (nrow(res.d7.Second_AM)*0.05), 
            ((table(res.d7.Second_AM$padj < 0.1)[2]) + (nrow(res.d7.Second_AM)*0.05)),
            ((table(res.d7.Second_AM$padj < 0.05)[2]) + (nrow(res.d7.Second_AM)*0.05)) ),
                  col=c("blue", "red", "red"), lty=c(1,2), lwd=c(1, 3)) # add line at 
dev.off()

# Write results csv
# group results 
resdata.d7.second.AM  <- merge(as.data.frame(res.d7.Second_AM_ordered), as.data.frame(counts(dds.d7, normalized=TRUE)), by="row.names", sort=FALSE) ## Merge with normalized count data
names(resdata.d7.second.AM)[1] <- "Gene"
resdata.d7.second.AM <- resdata.d7.second.AM[order(resdata.d7.second.AM$padj), ] # Order by adjusted p-value
path_out.d7 = 'C:/Users/samjg/Documents/My_Projects/Pgenerosa_TagSeq_Metabolomics/TagSeq/Analysis/Output/DESeq2/10cpm/Day7/' # call path
write.csv(resdata.d7.second.AM, paste(path_out.d7,"Day7.SecondDEGs_AvM_DESeq2results.csv")) # write

# volcano plot 
png("../Output/DESeq2/10cpm/Day7/Day7.Second.AM-VolcanoPlot.png", 1000, 1000, pointsize=20)
EnhancedVolcano(res.d7.Second_AM,
                lab = rownames(res.d7.Second_AM),
                x = 'log2FoldChange',
                y = 'padj',
                title = 'Day 7: Second Treatment (A v M)',
                subtitle = "DESeq2 - Differential expression",
                FCcutoff = 1,
                pCutoff = 0.05,
                pointSize = 4.0,
                labSize = 6.0,
                colAlpha = 1,
                legendPosition = 'right',
                col = c('black', '#E69F00', '#D55E00', '#0072B2'), # colorblindness pallette
                legendLabSize = 12,
                legendIconSize = 4.0,
                widthConnectors = 0.75)
dev.off()









# Secondary treatment A v S ========================================================================================================= 

# DEGs total for the primary treatment effect
d7.numDEGs_padj_S.AS               <- data.frame(table(res.d7.Second_AS$padj<0.05))[2,2] # DEGs - NOT considering LFC - just p adj
res.d7.Second_AS_ordered           <- res.d7.Second_AS[order(res.d7.Second_AS$padj), ] # Order by adjusted p-value - 110 DEGs
d7.num.UpReg_AS                    <- sum((res.d7.Second_AS_ordered$log2FoldChange[1:d7.numDEGs_padj_S.AS] > 0) == TRUE) #  LFC >= 1
d7.num.DownReg_AS                  <- sum((res.d7.Second_AS_ordered$log2FoldChange[1:d7.numDEGs_padj_S.AS] < 0) == TRUE) # LFC >= 1
d7.total_Sec.AS                    <- sum(d7.num.UpReg_AS,d7.num.DownReg_AS) # 106 TOTAL DEGs with the criteria pdj < 0.05 + LFC>1 (< -1)

resdata.d7.second.AS <- merge(as.data.frame(res.d7.Second_AS_ordered), as.data.frame(counts(dds.d7, normalized=TRUE)), by="row.names", sort=FALSE) ## Merge with normalized count data
names(resdata.d7.second.AS)[1] <- "Gene"
resdata.d7.second.AS <- resdata.d7.second.AS[order(resdata.d7.second.AS$padj), ] # Order by adjusted p-value
path_out.d7 = 'C:/Users/samjg/Documents/My_Projects/Pgenerosa_TagSeq_Metabolomics/TagSeq/Analysis/Output/DESeq2/10cpm/Day7/' # call path
write.csv(resdata.d7.second.AS, paste(path_out.d7,"Day7.SecondDEGs_AvS_DESeq2results.csv")) # write








# Secondary treatment M v S ========================================================================================================= 

# DEGs total for the primary treatment effect
d7.numDEGs_padj_MS               <- data.frame(table(res.d7.Second_MS$padj<0.05))[2,2] # DEGs - NOT considering LFC - just p adj
res.d7.Second_MS_ordered         <- res.d7.Second_MS[order(res.d7.Second_MS$padj), ] # Order by adjusted p-value - 110 DEGs
d7.num.UpReg_MS                  <- sum((res.d7.Second_MS_ordered$log2FoldChange[1:d7.numDEGs_padj_MS] > 0) == TRUE) #  LFC >= 1
d7.num.DownReg_MS                <- sum((res.d7.Second_MS_ordered$log2FoldChange[1:d7.numDEGs_padj_MS] < 0) == TRUE) # LFC >= 1
d7.total_Sec.MS                  <- sum(d7.num.UpReg_MS,d7.num.DownReg_MS) # 106 TOTAL DEGs with the criteria pdj < 0.05 + LFC>1 (< -1)

resdata.d7.second.MS <- merge(as.data.frame(res.d7.Second_MS_ordered), as.data.frame(counts(dds.d7, normalized=TRUE)), by="row.names", sort=FALSE) ## Merge with normalized count data
names(resdata.d7.second.MS)[1] <- "Gene"
resdata.d7.second.MS <- resdata.d7.second.MS[order(resdata.d7.second.MS$padj), ] # Order by adjusted p-value
path_out.d7 = 'C:/Users/samjg/Documents/My_Projects/Pgenerosa_TagSeq_Metabolomics/TagSeq/Analysis/Output/DESeq2/10cpm/Day7/' # call path
write.csv(resdata.d7.second.MS, paste(path_out.d7,"Day7.SecondDEGs_MvS_DESeq2results.csv")) # write



# primary x Secondary treatment  MA vs. AM  ============================================================================================ 
# view the comments above as to why to oouput this pairwise dataset specifically 
# in summary... the for loop later in this script outputs EACH full pariwise DE test and AM v. MA appears to carry the weight 
# of this UNIQUE treatment effect of 106 DEGs! (all other second and third treatment effects are minimal with 0 - 13 DEGs total)

# DEGs total for the primary treatment effect
d7.numDEGs_padj_MA_AM            <- data.frame(table(res.d7.Second_MA_AM$padj<0.05))[2,2] # DEGs - NOT considering LFC - just p adj
res.d7.MA_AM_ordered             <- res.d7.Second_MA_AM[order(res.d7.Second_MA_AM$padj), ] # Order by adjusted p-value - 110 DEGs
d7.num.UpReg_MA_AM               <- sum((res.d7.MA_AM_ordered$log2FoldChange[1:d7.numDEGs_padj_MA_AM] > 0) == TRUE) #  LFC >= 1
d7.num.DownReg_MA_AM             <- sum((res.d7.MA_AM_ordered$log2FoldChange[1:d7.numDEGs_padj_MA_AM] < 0) == TRUE) # LFC >= 1
d7.total_MA_AM                   <- sum(d7.num.UpReg_MA_AM,d7.num.DownReg_MA_AM) # 106 TOTAL DEGs with the criteria pdj < 0.05 + LFC>1 (< -1)

resdata.d7.MA_AM <- merge(as.data.frame(res.d7.MA_AM_ordered), as.data.frame(counts(dds.d7, normalized=TRUE)), by="row.names", sort=FALSE) ## Merge with normalized count data
names(resdata.d7.MA_AM)[1] <- "Gene"
resdata.d7.MA_AM <- resdata.d7.MA_AM[order(resdata.d7.MA_AM$padj), ] # Order by adjusted p-value
path_out.d7 = 'C:/Users/samjg/Documents/My_Projects/Pgenerosa_TagSeq_Metabolomics/TagSeq/Analysis/Output/DESeq2/10cpm/Day7/' # call path
write.csv(resdata.d7.MA_AM, paste(path_out.d7,"Day7.SecondDEGs_MA_AM_DESeq2results.csv")) # write




#MASTER MAIN EFFECTS TABLE FOR OUTPUT  ========================================================================================================
# rbind all significant DESeq results in main effect models 
Day7_MainEffectMaster <- rbind( 
  (as.data.frame(resdata.d7.primary)[c(1,3,7)]   %>% dplyr::mutate(Day_treatment = "Day7", Treatment = "Primary", Pairwise = "AvM") %>%  dplyr::filter(padj < 0.05) ),
  (as.data.frame(resdata.d7.second.AM)[c(1,3,7)] %>% dplyr::mutate(Day_treatment = "Day7", Treatment = "Second",  Pairwise = "AvM") %>%  dplyr::filter(padj < 0.05) ),
  (as.data.frame(resdata.d7.second.AS)[c(1,3,7)] %>% dplyr::mutate(Day_treatment = "Day7", Treatment = "Second",  Pairwise = "AvS") %>%  dplyr::filter(padj < 0.05) ),
  (as.data.frame(resdata.d7.second.MS)[c(1,3,7)] %>% dplyr::mutate(Day_treatment = "Day7", Treatment = "Second",  Pairwise = "MvS") %>%  dplyr::filter(padj < 0.05) ) )
colnames(Day7_MainEffectMaster)
# load annotation to merge by GeneID 
Pgen_annot <- read.delim(file="../../Seq_details/Panopea-generosa-genes-annotations.txt", header=F) 
Pgen_annot <- Pgen_annot[,c(1,5,7)]
names(Pgen_annot) <- c('Gene', 'Uniprot', 'Gene_terms_EC')
# merge the master with the Pgen_annot and save the csv file 
Day7_MainEffectMaster_2 <- merge(Day7_MainEffectMaster, Pgen_annot, by ="Gene")
path_out.d7 = 'C:/Users/samjg/Documents/My_Projects/Pgenerosa_TagSeq_Metabolomics/TagSeq/Analysis/Output/DESeq2/10cpm/Day7/' # call path to save
write.csv(Day7_MainEffectMaster_2, paste(path_out.d7,"Day7.AllMainEffectsDEGs_DESeq2results.csv")) # write








#  Data transformations for heatmap and PCA visuals  =======================================================================================


# Plot dispersions 
png("../Output/DESeq2/10cpm/Day7/Day7.dispersions.png", 1000, 1000, pointsize=20)
plotDispEsts(dds.d7, main="Day 7 dispersions")
dev.off()

# rlog - regularized log transformation of origin count data to log2 scale - fit for each sample and dist. of coefficients in the data
rlog.d7<- rlogTransformation(dds.d7) # rlog transform (regularized log)

png("../Output/DESeq2/10cpm/Day7/Day7.rlog_histogram.png", 1000, 1000, pointsize=20)# diagnostics of transformation # Histogram and sd plot
hist(assay(rlog.d7)) # view histogram 
dev.off()
png("../Output/DESeq2/10cpm/Day7/Day7.rlog_mean_sd.png", 1000, 1000, pointsize=20)
meanSdPlot(assay(rlog.d7)) # shows the sd y axis (sq root of varaince in all samples) - flat curve may seem like a goals, BUT may be unreasonable in cases with MANY true DEGs from experimental conditions
dev.off()



# PCA plot rlog ----------------------------------------------- #

pcaData_d7 <- plotPCA(rlog.d7, intgroup = c( "Primary_Treatment", "Second_Treament"), returnData = TRUE)
percentVar <- round(100 * attr(pcaData_d7, "percentVar"))
png("../Output/DESeq2/10cpm/Day7/Day7.rlog_PCA.png", 1000, 1000, pointsize=20)
PCA_plot_d14 <- ggplot(pcaData_d7, aes(x = PC1, y = PC2, color = Primary_Treatment, shape = Second_Treament)) +
                    scale_shape_manual(values = c(5, 19, 15)) +
                    geom_text(aes(label=name),hjust=0.2, vjust=1.4, size=3) +
                    geom_point(size =6) +
                    theme_minimal() +
                    theme(text = element_text(size=15)) +
                    ggtitle("PCA: Day7 (rlog)") +
                    xlab(paste0("PC1: ", percentVar[1], "% variance")) +
                    ylab(paste0("PC2: ", percentVar[2], "% variance")) +
                    coord_fixed()
PCA_plot_d14 + stat_ellipse(geom="polygon", aes(fill = Second_Treament), 
                 alpha = 0.2, 
                 show.legend = FALSE, 
                 level = 0.95) +  
                 theme(panel.grid = element_blank(), 
           panel.border = element_rect(fill= "transparent"))

dev.off()



# pheatmap rlog ----------------------------------------------- #



# COLORBLINDESS WHEEL
# Run pheatmap
# colorblind palette
# colorBlindBlack8  <- c("#000000", "#E69F00", "#56B4E9", "#009E73", 
#                        "#F0E442", "#0072B2", "#D55E00", "#CC79A7")
# pie(rep(1, 8), col = colorBlindBlack8)

save_pheatmap <- function(x, filename, width=1000, height=960) { # template for saving pheatmap outputs
  stopifnot(!missing(x))
  stopifnot(!missing(filename))
  png(filename,width = width, height=height)
  grid::grid.newpage()
  grid::grid.draw(x$gtable)
  dev.off()
}

# Primary Treatmetn M v A
topgenes.IDs <- head(rownames(res.d7_Prim_AvM_group.ordered),40) # call the first 40; FDR < 0.05 correct pdj
topgenes_rlog.counts <- assay(rlog.d7)[topgenes.IDs,]
topgenes_corrected <- topgenes_rlog.counts - rowMeans(topgenes_rlog.counts) # substract from the row mean to get +/- 0 to normalize a
df_annot.col <- as.data.frame(colData(dds.d7)[c("Primary_Treatment","Second_Treament")])

ann_colors = list(
     Primary_Treatment = c(A="#56B4E9", M="#E69F00"),
     Second_Treament = c(A = "#56B4E9", M = "#E69F00", S ="#D55E00"))

d7.rlog.heatmap <- pheatmap(topgenes_corrected, 
                            annotation_col=df_annot.col,
                            main = "Day 7 PRIMARY Ambient versus Moderate (106 total DEGs; 40 shown)",
                            #color = colorRampPalette(c("limegreen", "white", "gray39"))(50),
                            cutree_cols = 3,
                            cutree_rows = 2,
                            annotation_legend = TRUE,
                            annotation_colors = ann_colors,
                            show_rownames = FALSE,
                            labels_col=(paste(df_annot.col$Primary_Treatment,df_annot.col$Second_Treament, sep ="_")),
                            angle_col = 45,
                            fontsize = 8,
                            legend = TRUE)
#Save ppheatmap
save_pheatmap(d7.rlog.heatmap, filename = "../Output/DESeq2/10cpm/Day7/Day7.PRIMARY.rlog_heatmap.png")


# Secondary treatment A v M
topgenes.IDs <- head(rownames(res.d7.Second_AM_ordered),40) # call the first 40; FDR < 0.05 correct pdj
topgenes_rlog.counts <- assay(rlog.d7)[topgenes.IDs,]
topgenes_corrected <- topgenes_rlog.counts - rowMeans(topgenes_rlog.counts) # substract from the row mean to get +/- 0 to normalize a
df_annot.col <- as.data.frame(colData(dds.d7)[c("Primary_Treatment","Second_Treament")])

ann_colors = list(
     Primary_Treatment = c(A="#56B4E9", M="#E69F00"),
     Second_Treament = c(A = "#56B4E9", M = "#E69F00", S ="#D55E00"))

d7.rlog.heatmap.SEC <- pheatmap(topgenes_corrected, 
                            annotation_col=df_annot.col,
                            main = "Day 7 SECOND Ambient versus Moderate (106 total DEGs; 40 shown)",
                            #color = colorRampPalette(c("limegreen", "white", "gray39"))(50),
                            cutree_cols = 2,
                            cutree_rows = 2,
                            annotation_legend = TRUE,
                            annotation_colors = ann_colors,
                            show_rownames = FALSE,
                            labels_col=(paste(df_annot.col$Primary_Treatment,df_annot.col$Second_Treament, sep ="_")),
                            angle_col = 45,
                            fontsize = 8,
                            legend = TRUE)
#Save ppheatmap
save_pheatmap(d7.rlog.heatmap.SEC, filename = "../Output/DESeq2/10cpm/Day7/Day7.SECOND.AM.rlog_heatmap.png")


```

#Day 14
```{r Day14_Main_effects}
# TEST! (read below)
# 'Interaction' term tests - run this to see that the 'Primary Treatment' main effect in NOT represented in the full effect model
test1 <- results(dds.d14, contrast = c("Primary_Treatment", "M", "A"),  alpha= 0.05) # D/N addresses Primary M vs A in reference level Secondary!!! 
# thus the model above measures the pairwsie comparisons of MA vs. AA
table(test1$padj<0.05) #  5 DEGs
test2 <- results(dds.d14.group, contrast = list(("All_TreatmentMA"),("All_TreatmentAA")),  alpha= 0.05) # look here at the group - same genes are DEGs
table(test2$padj<0.05) #  5 DEGs (same as above)



# Primary Treatmetn M v A ========================================================================================================== #



# group by contrast
res.d14_primary <- results(dds.d14.group, contrast = list(c("All_TreatmentAA", "All_TreatmentAM", "All_TreatmentAS"),  c("All_TreatmentMA", "All_TreatmentMM", "All_TreatmentMS")),  alpha= 0.05)
table(res.d14_primary$padj<0.05) # 429 DEGs

# main effect model (mo interaction term)
res.d14.primary_main <- results(dds.d14.main, contrast = c("Primary_Treatment", "A", "M"),  alpha= 0.05)
table(res.d14.primary_main$padj<0.05) #  431 DEGs 

# full model 
res.d14.primary_fullmod <- results(dds.d14, contrast = c("Primary_Treatment", "A", "M"),  alpha= 0.05)
table(res.d14.primary_fullmod$padj<0.05) #  5 DEGs 

# Group contrast - this tests the full model effect above - appears the interaction term in fullmod is the pariwise effect under the reference level in treatment 2
res.d14.primary_grouptest<- results(dds.d14.group, contrast = list(("All_TreatmentAA"),("All_TreatmentMA")),  alpha= 0.05) # look here at the group - same genes are DEGs
table(res.d14.primary_grouptest$padj<0.05) #  5 DEGs (same as above)



# Secondary treatment A v M ========================================================================================================= #


# group by contrast
res.d14.Second_AM <- results(dds.d14.group, contrast = list(c("All_TreatmentAA", "All_TreatmentMA"),c("All_TreatmentAM", "All_TreatmentMM")),  alpha= 0.05)
table(res.d14.Second_AM$padj<0.05) # 0 DEGs

# main effect model (mo interaction term)
res.d14.Second_AM_main <- results(dds.d14.main, contrast = c("Second_Treament", "A", "M"),  alpha= 0.05)
table(res.d14.Second_AM_main$padj<0.05) #  0 DEGs 

# full model 
res.d14.Second_AM_fullmod <- results(dds.d14, contrast = c("Second_Treament", "A", "M"),  alpha= 0.05) # appears this is actually the pairwise comparisons between AA and AM
table(res.d14.Second_AM_fullmod$padj<0.05) #  0 DEGs 

# Group contrast - test the putative pariwise comparison with the ref level in fullmod
res.d14.Second_AM_grouptest<- results(dds.d14.group, contrast = list(("All_TreatmentAA"),("All_TreatmentAM")),  alpha= 0.05) # look here at the group constrasts
table(res.d14.Second_AM_grouptest$padj<0.05) #  0 DEGs 



# Secondary treatment A v S ========================================================================================================= #


# group by contrast
res.d14.Second_AS <- results(dds.d14.group, contrast = list(c("All_TreatmentAA", "All_TreatmentMA"),c("All_TreatmentAS", "All_TreatmentMS")),  alpha= 0.05)
table(res.d14.Second_AS$padj<0.05) #  0 DEGs

# main effect model (mo interaction term)
res.d14.Second_AS_main <- results(dds.d14.main, contrast = c("Second_Treament", "A", "S"),  alpha= 0.05)
table(res.d14.Second_AS_main$padj<0.05) #  0 DEGs 

# full model 
res.d14.Second_AS_fullmod <- results(dds.d14, contrast = c("Second_Treament", "A", "S"),  alpha= 0.05) # appears this is actually the pairwise comparisons between AA and AS
table(res.d14.Second_AS_fullmod$padj<0.05) #  1 DEGs 

# Group contrast - test the putative pariwise comparison with the ref level in fullmod
res.d14.Second_AS_grouptest<- results(dds.d14.group, contrast = list(("All_TreatmentAA"),("All_TreatmentAS")),  alpha= 0.05) # look here at the group constrasts
table(res.d14.Second_AS_grouptest$padj<0.05) #  1 DEGs (same as above)



# Secondary treatment M v S ========================================================================================================= #


# group by contrast
res.d14.Second_MS <- results(dds.d14.group, contrast = list(c("All_TreatmentAM", "All_TreatmentMM"),c("All_TreatmentAS", "All_TreatmentMS")),  alpha= 0.05)
table(res.d14.Second_MS$padj<0.05) #  0 DEGs

# main effect model (mo interaction term)
res.d14.Second_MS_main <- results(dds.d14.main, contrast = c("Second_Treament", "M", "S"),  alpha= 0.05)
table(res.d14.Second_MS_main$padj<0.05) #  0 DEGs 

# full model 
res.d14.Second_MS_fullmod <- results(dds.d14, contrast = c("Second_Treament", "M", "S"),  alpha= 0.05) # appears this is actually the pairwise comparisons between AM and AS
table(res.d14.Second_MS_fullmod$padj<0.05) #  0 DEGs 

# Group contrast - test the putative pariwise comparison with the ref level in fullmod
res.d14.Second_MS_grouptest<- results(dds.d14.group, contrast = list(("All_TreatmentAM"),("All_TreatmentAS")),  alpha= 0.05) # look here at the group constrasts
table(res.d14.Second_MS_grouptest$padj<0.05) #  0 DEGs 




# TABLES  &  Volcano Plots




# Primary Treatmetn M v A ========================================================================================================= #

# table(res.d14_P.AvM_group$padj<0.05) 
d14.numDEGs_padj               <-data.frame(table(res.d14_primary$padj<0.05))[2,2] # DEGs - NOT considering LFC - just p adj
res.d14_primary.ordered        <- res.d14_primary[order(res.d14_primary$padj), ] # Order by adjusted p-value
d14.num.UpReg                  <- sum((res.d14_primary.ordered$log2FoldChange[1:d14.numDEGs_padj] > 0) == TRUE) #  LFC >= 1
d14.num.DownReg                <- sum((res.d14_primary.ordered$log2FoldChange[1:d14.numDEGs_padj] < 0) == TRUE) # LFC >= 1
d14.total                      <- sum(d14.num.DownReg,d14.num.UpReg) # 429 sum of DEGs with the criteria pdj < 0.05 + LFC>1 (< -1) - 507

# Histogram
hist(res.d14_primary$pvalue, breaks=20, col="grey") # view histogram, 
abline(h=c( (nrow(res.d14_primary)*0.1), 
            ((table(res.d14_primary$padj < 0.1)[2]) + (nrow(res.d14_primary)*0.1)),
            ((table(res.d14_primary$padj < 0.05)[2]) + (nrow(res.d14_primary)*0.1)) ),
                  col=c("blue", "red", "red"), lty=c(1,2), lwd=c(1, 3)) # add line at 

# Write csv
# group model 
resdata.d14.primary  <- merge(as.data.frame(res.d14_primary.ordered), as.data.frame(counts(dds.d14.group, normalized=TRUE)), by="row.names", sort=FALSE) # Merge with normalized count
names(resdata.d14.primary)[1] <- "Gene"
resdata.d14.primary <- resdata.d14.primary[order(resdata.d14.primary$padj), ] # Order by adjusted p-value
path_out.d14 = 'C:/Users/samjg/Documents/My_Projects/Pgenerosa_TagSeq_Metabolomics/TagSeq/Analysis/Output/DESeq2/10cpm/Day14/' # call path
write.csv(resdata.d14.primary, paste(path_out.d14,"Day14.Primary_DESeq2results.csv")) # write

# main effect model
resdata.d14.primary_main  <- merge(as.data.frame(res.d14.primary_main), as.data.frame(counts(dds.d14, normalized=TRUE)), by="row.names", sort=FALSE) # Merge with normalized count
names(resdata.d14.primary_main)[1] <- "Gene"
resdata.d14.primary_main <- resdata.d14.primary_main[order(resdata.d14.primary_main$padj), ] # Order by adjusted p-value
path_out.d14.tests = 'C:/Users/samjg/Documents/My_Projects/Pgenerosa_TagSeq_Metabolomics/TagSeq/Analysis/Output/DESeq2/10cpm/Day14/D14_res_model_tests/' # call path
write.csv(resdata.d14.primary_main, paste(path_out.d14.tests,"Day14.Primary_res_main.mod.csv")) # write

# full model
resdata.d14.primary_fullmod  <- merge(as.data.frame(res.d14.primary_fullmod), as.data.frame(counts(dds.d14, normalized=TRUE)), by="row.names", sort=FALSE) # Merge with normalized count 
names(resdata.d14.primary_fullmod)[1] <- "Gene"
resdata.d14.primary_fullmod <- resdata.d14.primary_fullmod[order(resdata.d14.primary_fullmod$padj), ] # Order by adjusted p-value
path_out.d14.tests = 'C:/Users/samjg/Documents/My_Projects/Pgenerosa_TagSeq_Metabolomics/TagSeq/Analysis/Output/DESeq2/10cpm/Day14/D14_res_model_tests/' # call path
write.csv(resdata.d14.primary_fullmod, paste(path_out.d14.tests,"Day14.Primary_res_full.mod.csv")) # write

# group test of full model
resdata.d14.primary_grouptest  <- merge(as.data.frame(res.d14.primary_grouptest), as.data.frame(counts(dds.d14, normalized=TRUE)), by="row.names", sort=FALSE)# Merge with normalized 
names(resdata.d14.primary_grouptest)[1] <- "Gene"
resdata.d14.primary_grouptest <- resdata.d14.primary_grouptest[order(resdata.d14.primary_grouptest$padj), ] # Order by adjusted p-value
path_out.d14.tests = 'C:/Users/samjg/Documents/My_Projects/Pgenerosa_TagSeq_Metabolomics/TagSeq/Analysis/Output/DESeq2/10cpm/Day14/D14_res_model_tests/' # call path
write.csv(resdata.d14.primary_grouptest, paste(path_out.d14.tests,"Day14.Primary_res_grouptest.mod.csv")) # write

# volcano plot 
png("../Output/DESeq2/10cpm/Day14/Day14.Primary-VolcanoPlot.png", 1000, 1000, pointsize=20)
# Day14.volcano.primary <- 
EnhancedVolcano(res.d14_primary,
          lab = row.names(res.d14_primary),
          x = 'log2FoldChange',
          y = 'padj',
          title = 'Day 14: Primary Treatment (A v M)',
          subtitle = "DESeq2 - Differential expression",
          selectLab = c("PGEN_.00g057310", "PGEN_.00g279590", "PGEN_.00g092270", "PGEN_.00g092270", "PGEN_.00g266130"), # call persistant DEGs (day 7,14,21) with mean LFC > 5
          # selectLab = c("PGEN_.00g003940", "PGEN_.00g057310", "PGEN_.00g092270" ,"PGEN_.00g169480",
          #               "PGEN_.00g188250", "PGEN_.00g199810", "PGEN_.00g211250" ,"PGEN_.00g266130",
          #               "PGEN_.00g272710" ,"PGEN_.00g279590","PGEN_.00g348220", "PGEN_.00g348570"), # the list of genes from 'Const.DEGs.list'
                              xlab = bquote(~Log[2]~ 'fold change'),
                              ylab = bquote(~-Log[10]~ 'Padj'),
                              pCutoff = 0.05,
                              FCcutoff = 1.0,
                              pointSize = 3.0,
                              labSize = 3.0,
                              labCol = 'black',
                              labFace = 'bold',
                              boxedLabels = TRUE,
                            # parseLabels = TRUE,
                              col = c('black', '#E69F00', '#D55E00', '#0072B2'), # colorblindness pallette
                              colAlpha = 1,
                              legendPosition = 'right',
                              legendLabSize = 10,
                              legendIconSize = 3.0,
                              drawConnectors = TRUE,
                              widthConnectors = 0.2,
                              colConnectors = 'black')
dev.off()





#MASTER MAIN EFFECTS TABLE FOR OUTPUT  ========================================================================================================
# rbind all significant DESeq results in main effect models NOTE: only DEGs from the primary effect model 
Day14_MainEffectMaster <-  #rbind( 
  (as.data.frame(resdata.d14.primary)[c(1,3,7)]   %>% dplyr::mutate(Day_treatment = "Day14", Treatment = "Primary", Pairwise = "AvM") %>%  dplyr::filter(padj < 0.05) ) #,
  #(as.data.frame(resdata.d7.second.AM)[c(1,3,7)] %>% dplyr::mutate(Day_treatment = "Day7", Treatment = "Second",  Pairwise = "AvM") %>%  dplyr::filter(padj < 0.05) ),
  #(as.data.frame(resdata.d7.second.AS)[c(1,3,7)] %>% dplyr::mutate(Day_treatment = "Day7", Treatment = "Second",  Pairwise = "AvS") %>%  dplyr::filter(padj < 0.05) ),
  #(as.data.frame(resdata.d7.second.MS)[c(1,3,7)] %>% dplyr::mutate(Day_treatment = "Day7", Treatment = "Second",  Pairwise = "MvS") %>%  dplyr::filter(padj < 0.05) ) )
# load annotation to merge by GeneID 
Pgen_annot <- read.delim(file="../../Seq_details/Panopea-generosa-genes-annotations.txt", header=F) 
Pgen_annot <- Pgen_annot[,c(1,5,7)]
names(Pgen_annot) <- c('Gene', 'Uniprot', 'Gene_terms_EC')
# merge the master with the Pgen_annot and save the csv file 
Day14_MainEffectMaster_2 <- merge(Day14_MainEffectMaster, Pgen_annot, by ="Gene")
path_out.d14 = 'C:/Users/samjg/Documents/My_Projects/Pgenerosa_TagSeq_Metabolomics/TagSeq/Analysis/Output/DESeq2/10cpm/Day14/' # call path to save
write.csv(Day14_MainEffectMaster_2, paste(path_out.d14,"Day14.AllMainEffectsDEGs_DESeq2results.csv")) # write











#  Data transformations for heatmap and PCA visuals



# Plot dispersions 
png("../Output/DESeq2/10cpm/Day14/Day14.dispersions.png", 1000, 1000, pointsize=20)
plotDispEsts(dds.d14, main="Day 14 dispersions")
dev.off()

# Data transformations for heatmap and PCA visuals ======================================================================================= #
# rlog - regularized log transformation of origin count data to log2 scale - fit for each sample and dist. of coefficients in the data
rlog.d14<- rlogTransformation(dds.d14) # rlog transform (regularized log)

png("../Output/DESeq2/10cpm/Day14/Day14.rlog_histogram.png", 1000, 1000, pointsize=20)# diagnostics of transformation # Histogram and sd plot
hist(assay(rlog.d14)) # view histogram 
dev.off()
png("../Output/DESeq2/10cpm/Day14/Day14.rlog_mean_sd.png", 1000, 1000, pointsize=20)
meanSdPlot(assay(rlog.d14)) # shows the sd y axis (sq root of varaince in all samples) - flat curve may seem like a goals, BUT may be unreasonable in cases with MANY true DEGs from experimental conditions
dev.off()


# PCA plot rlog ----------------------------------------------- #

pcaData_d14 <- plotPCA(rlog.d14, intgroup = c( "Primary_Treatment", "Second_Treament"), returnData = TRUE)
percentVar <- round(100 * attr(pcaData_d14, "percentVar"))
png("../Output/DESeq2/10cpm/Day14/Day14.rlog_PCA.png", 1000, 1000, pointsize=20)
ggplot(pcaData_d14, aes(x = PC1, y = PC2, color = Primary_Treatment, shape = Second_Treament)) +
  scale_shape_manual(values = c(5, 19, 15)) +
  geom_text(aes(label=name),hjust=0.2, vjust=1.4, size=3) +
  geom_point(size =6) +
  theme_classic() +
  theme(text = element_text(size=15)) +
  ggtitle("PCA: Day14 (rlog)") +
  xlab(paste0("PC1: ", percentVar[1], "% variance")) +
  ylab(paste0("PC2: ", percentVar[2], "% variance")) +
  coord_fixed()
dev.off()




# pheatmap rlog ----------------------------------------------- #




# Plot heat map rlog------------------------ #
save_pheatmap <- function(x, filename, width=1000, height=960) { # template for saving pheatmap outputs
  stopifnot(!missing(x))
  stopifnot(!missing(filename))
  png(filename,width = width, height=height)
  grid::grid.newpage()
  grid::grid.draw(x$gtable)
  dev.off()
}
# parameters for the pheatmap
topgenes.IDs <- head(rownames(res.d14_primary.ordered),40) # call the first 13 genes that are differentially expressed, FDR < 0.05 correct pdj
topgenes_rlog.counts <- assay(rlog.d14)[topgenes.IDs,] 
topgenes_corrected <- topgenes_rlog.counts - rowMeans(topgenes_rlog.counts) # substract from the row mean to get +/- 0 to normalizze and ease the asthetic
df_annot.col <- as.data.frame(colData(dds.d14)[c("Primary_Treatment","Second_Treament")])


# colorblind palette
# colorBlindBlack8  <- c("#000000", "#E69F00", "#56B4E9", "#009E73", 
#                        "#F0E442", "#0072B2", "#D55E00", "#CC79A7")
# pie(rep(1, 8), col = colorBlindBlack8)

ann_colors = list(
    Primary_Treatment = c(A="#56B4E9", M="#E69F00"),
    Second_Treament = c(A = "#56B4E9", M = "#E69F00", S ="#D55E00"))

d14.rlog.heatmap <- pheatmap(topgenes_corrected, 
                            annotation_col=df_annot.col, 
                            main = "Day 14 PRIMARY Ambient versus Moderate (421 total DEGs; 40 shown)",
                            #color = colorRampPalette(c("limegreen", "white", "gray39"))(50),
                            cutree_cols = 2,
                            cutree_rows = 2,
                            annotation_legend = TRUE,
                            annotation_colors = ann_colors,
                            show_rownames = FALSE,
                            labels_col=(paste(df_annot.col$Primary_Treatment,df_annot.col$Second_Treament, sep ="_")), 
                            angle_col = 45,
                            fontsize = 8,
                            legend = TRUE)
#Save ppheatmap
save_pheatmap(d14.rlog.heatmap, filename = "../Output/DESeq2/10cpm/Day14/Day14.PRIMARY.rlog_heatmap.png")
```

#Day 21
```{r Day21_Main_effects}
# TEST! (read below)
# 'Interaction' term tests - run this to see that the 'Primary Treatment' main effect in NOT represented in the full effect model
test1 <- results(dds.d21, contrast = c("Primary_Treatment", "M", "A"),  alpha= 0.05) # D/N addresses Primary M vs A in reference level Secondary!!! 
# thus the model above measures the pairwsie comparisons of MA vs. AA
table(test1$padj<0.05) #  73 DEGs
test2 <- results(dds.d21.group, contrast = list(("All_TreatmentMAA"),("All_TreatmentAAA")),  alpha= 0.05) # look here at the group - same genes are DEGs
table(test2$padj<0.05) #  68 DEGs (same as above)




# Primary Treatmetn M v A ========================================================================================================== #

res.d21_primary <- results(dds.d21.group, contrast = list(c("All_TreatmentAAA", "All_TreatmentAAM", "All_TreatmentAMA", "All_TreatmentAMM", "All_TreatmentASA","All_TreatmentASM"), c("All_TreatmentMAA", "All_TreatmentMAM", "All_TreatmentMMA", "All_TreatmentMMM", "All_TreatmentMSA", "All_TreatmentMSM")),  alpha= 0.1)
table(res.d21_primary$padj<0.05) # 155 DEGs

# main effect model (mo interaction term)
res.d21.primary_main <- results(dds.d21.main, contrast = c("Primary_Treatment", "A", "M"),  alpha= 0.05)
table(res.d21.primary_main$padj<0.05) #  165 DEGs 

# full model 
res.d21.primary_fullmod <- results(dds.d21, contrast = c("Primary_Treatment", "A", "M"),  alpha= 0.05)
table(res.d21.primary_fullmod$padj<0.05) #  73 DEGs 

# Group contrast - this tests the full model effect above - appears the interaction term in fullmod is the pariwise effect under the reference level in treatment 2
res.d21.primary_grouptest<- results(dds.d21.group, contrast = list(("All_TreatmentAAA"),("All_TreatmentMAA")),  alpha= 0.05) # look here at the group - same genes are DEGs
table(res.d21.primary_grouptest$padj<0.05) #  68 DEGs (same as above)







# Secondary treatment A v M ========================================================================================================= #

res.d21.Second_AM <- results(dds.d21.group, contrast = list(c("All_TreatmentAAA", "All_TreatmentAAM", "All_TreatmentMAA", "All_TreatmentMAM"), c("All_TreatmentAMA", "All_TreatmentAMM", "All_TreatmentMMA", "All_TreatmentMMM")),  alpha= 0.1)
table(res.d21.Second_AM$padj<0.05) # 1 DEGs

d21.numDEGs_padj_AM            <- data.frame(table(res.d21.Second_AM$padj<0.05))[2,2] # DEGs - NOT considering LFC - just p adj
res.d21_AM.ordered             <- res.d21.Second_AM[order(res.d21.Second_AM$padj), ] # Order by adjusted p-value
d21.num.UpReg_AM               <- sum((res.d21.Second_AM$log2FoldChange > 0)[1] == TRUE) #  LFC >= 1
d21.num.DownReg_AM             <- sum((res.d21.Second_AM$log2FoldChange < 0)[1] == TRUE) # LFC >= 1


# main effect model (mo interaction term)
res.d21.Second_AM_main <- results(dds.d21.main, contrast = c("Second_Treament", "A", "M"),  alpha= 0.05)
table(res.d21.Second_AM_main$padj<0.05) #  1 DEGs 

# full model 
res.d21.Second_AM_fullmod <- results(dds.d21, contrast = c("Second_Treament", "A", "M"),  alpha= 0.05) # appears this is actually the pairwise comparisons between AAA and AMA
table(res.d21.Second_AM_fullmod$padj<0.05) #  0 DEGs 

# Group contrast - test the putative pariwise comparison with the ref level in fullmod
res.d21.Second_AM_grouptest<- results(dds.d21.group, contrast = list(("All_TreatmentAAA"),("All_TreatmentAMA")),  alpha= 0.05) # look here at the group constrasts
table(res.d21.Second_AM_grouptest$padj<0.05) #  0 DEGs 








# Secondary treatment A v S ========================================================================================================= #

res.d21.Second_AS <- results(dds.d21.group, contrast = list(c("All_TreatmentAAA", "All_TreatmentAAM", "All_TreatmentMAA", "All_TreatmentMAM"), c("All_TreatmentASA", "All_TreatmentASM", "All_TreatmentMSA", "All_TreatmentMSM")),  alpha= 0.1)
table(res.d21.Second_AS$padj<0.05) # 1 DEGs

d21.numDEGs_padj_AS            <- data.frame(table(res.d21.Second_AS$padj<0.05))[2,2] # DEGs - NOT considering LFC - just p adj
res.d21_AS.ordered             <- res.d21.Second_AS[order(res.d21.Second_AS$padj), ] # Order by adjusted p-value
d21.num.UpReg_AS               <- sum((res.d21.Second_AS$log2FoldChange > 0)[1] == TRUE) #  LFC >= 1
d21.num.DownReg_AS             <- sum((res.d21.Second_AS$log2FoldChange < 0)[1] == TRUE) # LFC >= 1

# main effect model (mo interaction term)
res.d21.Second_AS_main <- results(dds.d21.main, contrast = c("Second_Treament", "A", "S"),  alpha= 0.05)
table(res.d21.Second_AS_main$padj<0.05) #  0 DEGs 

# full model 
res.d21.Second_AS_fullmod <- results(dds.d21, contrast = c("Second_Treament", "A", "S"),  alpha= 0.05) # appears this is actually the pairwise comparisons between AAA and ASA
table(res.d21.Second_AS_fullmod$padj<0.05) #  1 DEGs 

# Group contrast - test the putative pariwise comparison with the ref level in fullmod
res.d21.Second_AS_grouptest<- results(dds.d21.group, contrast = list(("All_TreatmentAAA"),("All_TreatmentASA")),  alpha= 0.05) # look here at the group constrasts
table(res.d21.Second_AS_grouptest$padj<0.05) #  1 DEGs (same as above)







# Secondary treatment M v S ========================================================================================================= #

res.d21.Second_MS <- results(dds.d21.group, contrast = list(c("All_TreatmentAMA", "All_TreatmentAMM", "All_TreatmentMMA", "All_TreatmentMMM"), c("All_TreatmentASA", "All_TreatmentASM", "All_TreatmentMSA", "All_TreatmentMSM")),  alpha= 0.1)
table(res.d21.Second_MS$padj<0.05) # 0 DEGs

# main effect model (mo interaction term)
res.d21.Second_MS_main <- results(dds.d21.main, contrast = c("Second_Treament", "M", "S"),  alpha= 0.05)
table(res.d21.Second_MS_main$padj<0.05) #  0 DEGs 

# full model 
res.d21.Second_MS_fullmod <- results(dds.d21, contrast = c("Second_Treament", "M", "S"),  alpha= 0.05) # appears this is actually the pairwise comparisons between AMA and ASA
table(res.d21.Second_MS_fullmod$padj<0.05) #  1 DEGs 

# Group contrast - test the putative pariwise comparison with the ref level in fullmod
res.d21.Second_MS_grouptest<- results(dds.d21.group, contrast = list(("All_TreatmentAMA"),("All_TreatmentASA")),  alpha= 0.05) # look here at the group constrasts
table(res.d21.Second_MS_grouptest$padj<0.05) #  2 DEGs 






# Third treatment A v M ========================================================================================================= #

res.d21.Third <- results(dds.d21.group, contrast = list(c("All_TreatmentAAA", "All_TreatmentAMA", "All_TreatmentASA", "All_TreatmentMAA", "All_TreatmentMMA", "All_TreatmentMSA"), c("All_TreatmentAAM", "All_TreatmentAMM", "All_TreatmentASM", "All_TreatmentMAM", "All_TreatmentMMM", "All_TreatmentMSM")),  alpha= 0.1)
table(res.d21.Third$padj<0.05) # 0 DEGs

# main effect model (mo interaction term)
res.d21.Third_main <- results(dds.d21.main, contrast = c("Third_Treatment", "A", "M"),  alpha= 0.05)
table(res.d21.Third_main$padj<0.05) #  0 DEGs 

# full model 
res.d21.Third_fullmod <- results(dds.d21, contrast = c("Third_Treatment", "A", "M"),  alpha= 0.05) # appears this is actually the pairwise comparisons between AAA and AAM
table(res.d21.Third_fullmod$padj<0.05) #  0 DEGs 

# Group contrast - test the putative pariwise comparison with the ref level in fullmod
res.d21.Third_grouptest<- results(dds.d21.group, contrast = list(("All_TreatmentAAA"),("All_TreatmentAAM")),  alpha= 0.05) # look here at the group constrasts
table(res.d21.Third_grouptest$padj<0.05) #  0 DEGs 






# TABLES  &  Volcano Plots - only DEGs found for primary AvM and Second AvS (1 DEG) and second AvM (1 DEG) 






# Primary Treatmetn M v A ========================================================================================================= #

d21.numDEGs_padj               <- data.frame(table(res.d21_primary$padj<0.05))[2,2] # DEGs - NOT considering LFC - just p adj
res.d21_primary.ordered        <- res.d21_primary[order(res.d21_primary$padj), ] # Order by adjusted p-value
d21.num.UpReg                  <- sum((res.d21_primary.ordered$log2FoldChange[1:d21.numDEGs_padj] > 0) == TRUE) #  LFC >= 1
d21.num.DownReg                <- sum((res.d21_primary.ordered$log2FoldChange[1:d21.numDEGs_padj] < 0) == TRUE) # LFC >= 1
d21.total                      <- sum(d21.num.DownReg,d21.num.UpReg) # 155 sum of DEGs with the criteria pdj < 0.05 + LFC>1 (< -1) - 230 genes
# Histogram
hist(res.d21_Prim_AvM_group$pvalue, breaks=20, col="grey") # view histogram,  
abline(h=c( (nrow(res.d21_Prim_AvM_group)*0.1), 
            ((table(res.d21_Prim_AvM_group$padj < 0.1)[2]) + (nrow(res.d21_Prim_AvM_group)*0.1)),
            ((table(res.d21_Prim_AvM_group$padj < 0.05)[2]) + (nrow(res.d21_Prim_AvM_group)*0.1)) ),
                  col=c("blue", "red", "red"), lty=c(1,2), lwd=c(1, 3)) # add line at 

# Write csv
resdata.d21.primary  <- merge(as.data.frame(res.d21_primary.ordered), as.data.frame(counts(dds.d21.group, normalized=TRUE)), by="row.names", sort=FALSE) ## Merge with normalized count data
names(resdata.d21.primary)[1] <- "Gene"
resdata.d21.primary <- resdata.d21.primary[order(resdata.d21.primary$padj), ] # Order by adjusted p-value
path_out.d21 = 'C:/Users/samjg/Documents/My_Projects/Pgenerosa_TagSeq_Metabolomics/TagSeq/Analysis/Output/DESeq2/10cpm/Day21/' # call path
write.csv(resdata.d21.primary, paste(path_out.d21,"Day21.Primary_DESeq2results.csv")) # write

# volcano plot 
png("../Output/DESeq2/10cpm/Day21/Day21.PrimaryTreatment-VolcanoPlot.png", 1000, 1000, pointsize=20)
EnhancedVolcano(res.d21_primary,
                lab = rownames(res.d21_primary),
                x = 'log2FoldChange',
                y = 'padj',
                title = 'Day 21: Primary Treatment (A v M)',
                subtitle = "DESeq2 - Differential expression",
                selectLab = c("PGEN_.00g057310", "PGEN_.00g279590", "PGEN_.00g092270", "PGEN_.00g092270", "PGEN_.00g266130"), # call persistant DEGs (day 7,14,21) with mean LFC > 5
                # selectLab = c("PGEN_.00g003940", "PGEN_.00g057310", "PGEN_.00g092270" ,"PGEN_.00g169480",
                #               "PGEN_.00g188250", "PGEN_.00g199810", "PGEN_.00g211250" ,"PGEN_.00g266130",
                #               "PGEN_.00g272710" ,"PGEN_.00g279590","PGEN_.00g348220", "PGEN_.00g348570"), # the list of genes from 'Const.DEGs.list'
                              xlab = bquote(~Log[2]~ 'fold change'),
                              ylab = bquote(~-Log[10]~ 'Padj'),
                              pCutoff = 0.05,
                              FCcutoff = 1.0,
                              pointSize = 3.0,
                              labSize = 3.0,
                              labCol = 'black',
                              labFace = 'bold',
                              boxedLabels = TRUE,
                            # parseLabels = TRUE,
                              col = c('black', '#E69F00', '#D55E00', '#0072B2'), # colorblindness pallette
                              colAlpha = 1,
                              legendPosition = 'right',
                              legendLabSize = 10,
                              legendIconSize = 3.0,
                              drawConnectors = TRUE,
                              widthConnectors = 0.2,
                              colConnectors = 'black')
dev.off()







# Second Treatmetn A v M  ========================================================================================================= #
resdata.d21.second_AM  <- merge(as.data.frame(res.d21_AM.ordered), as.data.frame(counts(dds.d21.group, normalized=TRUE)), by="row.names", sort=FALSE) #Merge with normalized count data
names(resdata.d21.second_AM)[1] <- "Gene"
resdata.d21.second_AM <- resdata.d21.second_AM[order(res.d21_AM.ordered$padj), ]






# Second Treatmetn A v S  ========================================================================================================= #
resdata.d21.second_AS  <- merge(as.data.frame(res.d21_AS.ordered), as.data.frame(counts(dds.d21.group, normalized=TRUE)), by="row.names", sort=FALSE) #Merge with normalized count data
names(resdata.d21.second_AS)[1] <- "Gene"
resdata.d21.second_AS <- resdata.d21.second_AS[order(resdata.d21.second_AS$padj), ]




#MASTER MAIN EFFECTS TABLE FOR OUTPUT  ========================================================================================================
# rbind all significant DESeq results in main effect models NOTE: only DEGs from the primary effect model 
Day21_MainEffectMaster <-  rbind( 
  (as.data.frame(resdata.d21.primary)[c(1,3,7)]   %>% dplyr::mutate(Day_treatment = "Day21", Treatment = "Primary", Pairwise = "AvM") %>%  dplyr::filter(padj < 0.05) ),
  (as.data.frame(resdata.d21.second_AM)[c(1,3,7)] %>% dplyr::mutate(Day_treatment = "Day21", Treatment = "Second",  Pairwise = "AvM") %>%  dplyr::filter(padj < 0.05) ),
  (as.data.frame(resdata.d21.second_AS)[c(1,3,7)] %>% dplyr::mutate(Day_treatment = "Day21", Treatment = "Second",  Pairwise = "AvS") %>%  dplyr::filter(padj < 0.05) ) )#,
  #(as.data.frame(resdata.d7.second.MS)[c(1,3,7)] %>% dplyr::mutate(Day_treatment = "Day7", Treatment = "Second",  Pairwise = "MvS") %>%  dplyr::filter(padj < 0.05) ) )
# load annotation to merge by GeneID 
Pgen_annot <- read.delim(file="../../Seq_details/Panopea-generosa-genes-annotations.txt", header=F) 
Pgen_annot <- Pgen_annot[,c(1,5,7)]
names(Pgen_annot) <- c('Gene', 'Uniprot', 'Gene_terms_EC')
# merge the master with the Pgen_annot and save the csv file 
Day21_MainEffectMaster_2 <- merge(Day21_MainEffectMaster, Pgen_annot, by ="Gene")
path_out.d21 = 'C:/Users/samjg/Documents/My_Projects/Pgenerosa_TagSeq_Metabolomics/TagSeq/Analysis/Output/DESeq2/10cpm/Day21/' # call path to save
write.csv(Day21_MainEffectMaster_2, paste(path_out.d21,"Day21.AllMainEffectsDEGs_DESeq2results.csv")) # write










#  Data transformations for heatmap and PCA visuals



# Plot dispersions 
png("../Output/DESeq2/10cpm/Day21/Day21.dispersions.png", 1000, 1000, pointsize=20)
plotDispEsts(dds.d21, main="Day 21 dispersions")
dev.off()

# Data transformations for heatmap and PCA visuals ======================================================================================= #
# rlog - regularized log transformation of origin count data to log2 scale - fit for each sample and dist. of coefficients in the data
rlog.d21<- rlogTransformation(dds.d21) # rlog transform (regularized log)

png("../Output/DESeq2/10cpm/Day21/Day21.rlog_histogram.png", 1000, 1000, pointsize=20)# diagnostics of transformation # Histogram and sd plot
hist(assay(rlog.d21)) # view histogram 
dev.off()
png("../Output/DESeq2/10cpm/Day21/Day21.rlog_mean_sd.png", 1000, 1000, pointsize=20)
meanSdPlot(assay(rlog.d21)) # shows the sd y axis (sq root of varaince in all samples) - flat curve may seem like a goals, BUT may be unreasonable in cases with MANY true DEGs from experimental conditions
dev.off()



# PCA plot rlog ----------------------------------------------- #

pcaData_d21 <- plotPCA(rlog.d21, intgroup = c( "Primary_Treatment", "Second_Treament", "Third_Treatment"), returnData = TRUE)
percentVar <- round(100 * attr(pcaData_d21, "percentVar"))
png("../Output/DESeq2/10cpm/Day21/Day21.rlog_PCA.png", 1000, 1000, pointsize=20)
ggplot(pcaData_d21, aes(x = PC1, y = PC2, color = Primary_Treatment, shape = Second_Treament, fill = Third_Treatment)) +
  scale_shape_manual(values = c(4, 19, 17)) +
  geom_text(aes(label=name),hjust=0.2, vjust=1.4, size=3) +
  geom_point(size =6) +
  theme_classic() +
  theme(text = element_text(size=15)) +
  ggtitle("PCA: Day21 (rlog)") +
  xlab(paste0("PC1: ", percentVar[1], "% variance")) +
  ylab(paste0("PC2: ", percentVar[2], "% variance")) +
  coord_fixed()
dev.off()



# pheatmap rlog ----------------------------------------------- #




# Plot heat map rlog------------------------ #
save_pheatmap <- function(x, filename, width=1000, height=960) { # template for saving pheatmap outputs
  stopifnot(!missing(x))
  stopifnot(!missing(filename))
  png(filename,width = width, height=height)
  grid::grid.newpage()
  grid::grid.draw(x$gtable)
  dev.off()
}
# parameters for the pheatmap
topgenes.IDs <- head(rownames(res.d21_primary.ordered),40) # call the first 13 genes that are differentially expressed, FDR < 0.05 correct pdj
topgenes_rlog.counts <- assay(rlog.d21)[topgenes.IDs,] 
topgenes_corrected <- topgenes_rlog.counts - rowMeans(topgenes_rlog.counts) # substract from the row mean to get +/- 0 to normalizze and ease the asthetic
df_annot.col <- as.data.frame(colData(dds.d21)[c("Primary_Treatment","Second_Treament","Third_Treatment")])

# colorblind palette
# colorBlindBlack8  <- c("#000000", "#E69F00", "#56B4E9", "#009E73", 
#                        "#F0E442", "#0072B2", "#D55E00", "#CC79A7")
# pie(rep(1, 8), col = colorBlindBlack8)

ann_colors = list(
    Primary_Treatment = c(A="#56B4E9", M="#E69F00"),
    Second_Treament = c(A = "#56B4E9", M = "#E69F00", S ="#D55E00"),
    Third_Treatment = c(A = "#56B4E9", M = "#E69F00"))

d21.rlog.heatmap <- pheatmap(topgenes_corrected, 
                            annotation_col=df_annot.col, 
                            main = "Day 21 PRIMARY Ambient versus Moderate (155 total DEGs; 40 shown)",
                            #color = colorRampPalette(c("limegreen", "white", "gray39"))(50),
                            cutree_cols = 3,
                            cutree_rows = 2,
                            annotation_legend = TRUE,
                            annotation_colors = ann_colors,
                            show_rownames = FALSE,
                            labels_col=(paste(df_annot.col$Primary_Treatment,df_annot.col$Second_Treament, df_annot.col$Third_Treatment)), 
                            angle_col = 45,
                            fontsize = 6,
                            legend = TRUE)
#Save ppheatmap
save_pheatmap(d21.rlog.heatmap, filename = "../Output/DESeq2/10cpm/Day21/Day21.PRIMARY.rlog_heatmap.png")
```

#Subseq Encounters - Pariwise Targets
```{r Subseqent Pairwise encounters}
# GROUP interactions - Biological Hypothesis testing ========================================================== #
# Hypothesis does previous history 'prime' response to subseqent stressor 
# Thus we will investigate AM v MM and MS v MS

# call path for all output cvs of DE results
path_out.subseq_intxns = 'C:/Users/samjg/Documents/My_Projects/Pgenerosa_TagSeq_Metabolomics/TagSeq/Analysis/Output/DESeq2/10cpm/DEGs_pairwise_subseqent_exposures/' # call path



# DAY 7



res.d7.AMvMM                  <- results(dds.d7.group, contrast = list(("All_TreatmentAM"),("All_TreatmentMM")),  alpha= 0.05)
table(res.d7.AMvMM$padj<0.05) #  5 DEGs
resdata.d7.AMvMM$geneSymbol   <- row.names(resdata.d7.AMvMM)
resdata.d7.AMvMM_Annot        <- merge(Pgen_annot, resdata.d7.AMvMM, by.x='V1',by.y='geneSymbol')
resdata.d7.AMvMM_Annot        <- as.data.frame(resdata.d7.AMvMM_Annot[order(resdata.d7.AMvMM_Annot$padj), ]) # Order by adjusted p-value
resdata.d7.AMvMM_Annot_DEGs   <- resdata.d7.AMvMM_Annot %>% filter(padj <0.05)
write.csv(resdata.d7.AMvMM_Annot_DEGs, paste(path_out.subseq_intxns,"Day7.AMvMM.csv")) # write

res.d7.ASvMS                 <- results(dds.d7.group, contrast = list(("All_TreatmentAS"),("All_TreatmentMS")),  alpha= 0.05)
table(res.d7.ASvMS$padj<0.05) #  15 DEGs
resdata.d7.ASvMS             <- as.data.frame(res.d7.ASvMS) # Order by adjusted p-value
resdata.d7.ASvMS$geneSymbol  <- row.names(resdata.d7.ASvMS)
res.d7.ASvMS_Annot           <- merge(Pgen_annot, resdata.d7.ASvMS, by.x='V1',by.y='geneSymbol')
res.d7.ASvMS_Annot           <- as.data.frame(res.d7.ASvMS_Annot[order(res.d7.ASvMS_Annot$padj), ]) # Order by adjusted p-value
res.d7.ASvMS_Annot_DEGs      <- res.d7.ASvMS_Annot %>% filter(padj <0.05)
write.csv(res.d7.ASvMS_Annot_DEGs, paste(path_out.subseq_intxns,"Day7.ASvMS.csv")) # write



# DAY 14



res.d14.AMvMM                  <- results(dds.d14.group, contrast = list(("All_TreatmentAM"),("All_TreatmentMM")),  alpha= 0.05)
table(res.d14.AMvMM$padj<0.05) #  5 DEGs
resdata.d14.AMvMM              <- as.data.frame(res.d14.AMvMM) # Order by adjusted p-value
resdata.d14.AMvMM$geneSymbol   <- row.names(resdata.d7.AMvMM)
resdata.d14.AMvMM_Annot        <- merge(Pgen_annot, resdata.d14.AMvMM, by.x='V1',by.y='geneSymbol')
resdata.d14.AMvMM_Annot        <- as.data.frame(resdata.d14.AMvMM_Annot[order(resdata.d14.AMvMM_Annot$padj), ]) # Order by adjusted p-value
resdata.d14.AMvMM_Annot_DEGs   <- resdata.d14.AMvMM_Annot %>% filter(padj <0.05)
write.csv(resdata.d14.AMvMM_Annot_DEGs, paste(path_out.subseq_intxns,"Day14.AMvMM.csv")) # write

res.d14.ASvMS            <- results(dds.d14.group, contrast = list(("All_TreatmentAS"),("All_TreatmentMS")),  alpha= 0.05)
table(res.d14.ASvMS$padj<0.05) #  15 DEGs
resdata.d14.ASvMS              <- as.data.frame(res.d14.ASvMS) # Order by adjusted p-value
resdata.d14.ASvMS             <- as.data.frame(resdata.d14.ASvMS) # Order by adjusted p-value
resdata.d14.ASvMS$geneSymbol  <- row.names(resdata.d14.ASvMS)
res.d14.ASvMS_Annot           <- merge(Pgen_annot, resdata.d14.ASvMS, by.x='V1',by.y='geneSymbol')
res.d14.ASvMS_Annot           <- as.data.frame(res.d14.ASvMS_Annot[order(res.d14.ASvMS_Annot$padj), ]) # Order by adjusted p-value
res.d14.ASvMS_Annot_DEGs      <- res.d14.ASvMS_Annot %>% filter(padj <0.05)
write.csv(res.d14.ASvMS_Annot_DEGs, paste(path_out.subseq_intxns,"Day14.ASvMS.csv")) # write




# DAY 21




res.d21.AMMvMMM                 <- results(dds.d21.group, contrast = list(("All_TreatmentAMM"),("All_TreatmentMMM")),  alpha= 0.05)
table(res.d21.AMMvMMM$padj<0.05) #   26 DEGs
resdata.d21.AMMvMMM              <- as.data.frame(res.d21.AMMvMMM) # Order by adjusted p-value
resdata.d21.AMMvMMM$geneSymbol   <- row.names(resdata.d21.AMMvMMM)
resdata.d21.AMMvMMM_Annot        <- merge(Pgen_annot, resdata.d21.AMMvMMM, by.x='V1',by.y='geneSymbol')
resdata.d21.AMMvMMM_Annot        <- as.data.frame(resdata.d21.AMMvMMM_Annot[order(resdata.d21.AMMvMMM_Annot$padj), ]) # Order by adjusted p-value
resdata.d21.AMMvMMM_Annot_DEGs   <- resdata.d21.AMMvMMM_Annot %>% filter(padj <0.05)
write.csv(resdata.d21.AMMvMMM_Annot_DEGs, paste(path_out.subseq_intxns,"Day21.AMMvMMM.csv")) # write

res.d21.ASMvMSM                <- results(dds.d21.group, contrast = list(("All_TreatmentASM"),("All_TreatmentMSM")),  alpha= 0.05)
table(res.d21.ASMvMSM$padj<0.05) #   26 DEGs
resdata.d21.ASMvMSM              <- as.data.frame(res.d21.ASMvMSM) # Order by adjusted p-value
resdata.d21.ASMvMSM$geneSymbol   <- row.names(resdata.d21.ASMvMSM)
resdata.d21.ASMvMSM_Annot        <- merge(Pgen_annot, resdata.d21.ASMvMSM, by.x='V1',by.y='geneSymbol')
resdata.d21.ASMvMSM_Annot        <- as.data.frame(resdata.d21.ASMvMSM_Annot[order(resdata.d21.ASMvMSM_Annot$padj), ]) # Order by adjusted p-value
resdata.d21.ASMvMSM_Annot_DEGs   <- resdata.d21.ASMvMSM_Annot %>% filter(padj <0.05)
write.csv(resdata.d21.ASMvMSM_Annot_DEGs, paste(path_out.subseq_intxns,"Day21.ASMvMSM.csv")) # write

```
# grid volcano - primary treatment effect with constant DEGs > 5 mean LFC
```{r grid volcano}
volcano.Grid.Primary <- grid.arrange(Day7.volcano.primary, Day14.volcano.primary,ncol=1, nrow=2, clip="off")

pdf("../Output/DESeq2/10cpm/volcano.Grid.Primary.pdf")
Day7.volcano.primary
# grid.arrange(Day7.volcano.primary, Day14.volcano.primary,ncol=1, nrow=2, clip="off")
dev.off()
```

# summary table of DEGs and occurances ("persistant" DEGs)
```{r Output_tables_and_knitr_summary_tables}
# write to table
# prep a table for knitr
Primary.DEGs.table <- data.frame(matrix(nrow = 4, ncol = 6)) # create a new data table
colnames(Primary.DEGs.table)<-c('Day', 'Treatment','Pairwise.compare', 'DEGs_total', 'DEGs_UP', 'DEGs_DOWN') 
Primary.DEGs.table$Day <- c("0", "7", "14", "21")
Primary.DEGs.table$Treatment <- "Primary"
Primary.DEGs.table$Pairwise.compare <- "A_v_M"
Primary.DEGs.table$DEGs_total <- c(d0.total, d7.total, d14.total, d21.total)
Primary.DEGs.table$DEGs_UP <- c(d0.num.UpReg, d7.num.UpReg, d14.num.UpReg, d21.num.UpReg)
Primary.DEGs.table$DEGs_DOWN <- c(d0.num.DownReg, d7.num.DownReg, d14.num.DownReg, d21.num.DownReg)


# FDR histograms - red line shows the false discovery rate of 0.05 called in each model as "alpha = 0.05"
# peaks at pdj < 0.05 above this line are considering true in the models
# hist(resd0.primary$pvalue, breaks=20, col="grey")          + abline(h=(nrow(resd0.primary)*0.05),col="red")
# hist(res.d7.primary$pvalue, breaks=20, col="grey")     + abline(h=(nrow(res.d7.primary)*0.05),col="red")
# hist(res.d14_P.AvM_group$pvalue, breaks=20, col="grey")    + abline(h=(nrow(res.d14_P.AvM_group)*0.05),col="red")
# hist(res.d21_Prim_AvM_group$pvalue, breaks=20, col="grey") + abline(h=(nrow(res.d21_Prim_AvM_group)*0.05),col="red")

# output tables for DEGs 
# Day 0 
DE_d0.primary <- merge(as.data.frame(resd0.primary), as.data.frame(counts(dds.d0, normalized=FALSE)), by="row.names", sort=FALSE) %>% dplyr::filter(padj < 0.05) %>% mutate(up = log2FoldChange < 0) %>% mutate(down = log2FoldChange > 0) # remember! this model is M vA whereas the rest ate A v M
write.csv(DE_d0.primary, "../Output/DESeq2/10cpm/Day0/DE_Day0_Primary.csv")

# Day 7 
DE_d7.primary <- merge(as.data.frame(res.d7.primary), as.data.frame(counts(dds.d7, normalized=FALSE)), by="row.names", sort=FALSE) %>% dplyr::filter(padj < 0.05) %>% mutate(up = log2FoldChange > 0) %>% mutate(down = log2FoldChange < 0)
write.csv(DE_d7.primary, "../Output/DESeq2/10cpm/Day7/DE_Day7_Primary.csv")

DE_d7.second <- merge(as.data.frame(res.d7.Second_AM), as.data.frame(counts(dds.d7, normalized=FALSE)), by="row.names", sort=FALSE) %>% dplyr::filter(padj < 0.05) %>% mutate(up = log2FoldChange > 0) %>% mutate(down = log2FoldChange < 0)
write.csv(DE_d7.second, "../Output/DESeq2/10cpm/Day7/DE_Day7_Second.csv") # second treatment A v M (106 DEGs)

DE_d7.MA_AM <- merge(as.data.frame(res.d7.Second_MA_AM), as.data.frame(counts(dds.d7, normalized=FALSE)), by="row.names", sort=FALSE) %>% dplyr::filter(padj < 0.05) %>% mutate(up = log2FoldChange > 0) %>% mutate(down = log2FoldChange < 0)
write.csv(DE_d7.MA_AM, "../Output/DESeq2/10cpm/Day7/DE_Day7_MA_AM.csv") # full model AM v MA (168 DEGs)


# Day 14 
DE_d14.primary <- merge(as.data.frame(res.d14_primary), as.data.frame(counts(dds.d14, normalized=FALSE)), by="row.names", sort=FALSE) %>% dplyr::filter(padj < 0.05) %>% mutate(up = log2FoldChange > 0) %>% mutate(down = log2FoldChange < 0)
write.csv(DE_d14.primary, "../Output/DESeq2/10cpm/Day14/DE_Day14_Primary.csv")

# Day 21
DE_d21.primary <- merge(as.data.frame(res.d21_primary), as.data.frame(counts(dds.d21, normalized=FALSE)), by="row.names", sort=FALSE) %>% dplyr::filter(padj < 0.05) %>% mutate(up = log2FoldChange > 0) %>% mutate(down = log2FoldChange < 0)
write.csv(DE_d21.primary, "../Output/DESeq2/10cpm/Day21/DE_Day21_Primary.csv")


d0_P.DEGs_boolean <- DE_d0.primary %>% dplyr::select(c("Row.names"
,"log2FoldChange","padj","up","down")) %>% dplyr::mutate(Day_Trmt = "Day0_Primary_AvM")

d7_P.DEGs_boolean <- DE_d7.primary %>% dplyr::select(c("Row.names"
,"log2FoldChange","padj","up","down")) %>% dplyr::mutate(Day_Trmt = "Day7_Primary_AvM")

d14_P.DEGs_boolean <- DE_d14.primary %>% dplyr::select(c("Row.names"
,"log2FoldChange","padj","up","down")) %>% dplyr::mutate(Day_Trmt = "Day14_Primary_AvM")

d21_P.DEGs_boolean <- DE_d21.primary %>% dplyr::select(c("Row.names"
,"log2FoldChange","padj","up","down")) %>% dplyr::mutate(Day_Trmt = "Day21_Primary_AvM")

All_DEGs_Primary.Ambieant_v_Moderate <- rbind(d0_P.DEGs_boolean,d7_P.DEGs_boolean, d14_P.DEGs_boolean,d21_P.DEGs_boolean)
write.csv(All_DEGs_Primary.Ambieant_v_Moderate, "../Output/DESeq2/10cpm/DE_PrimaryTreatment.All.csv")

# Summary table for the frequency that DEGs occured 
# unique(All_DEGs_Primary.Ambieant_v_Moderate$Day_Trmt) # unique DAy_treat to filter
DE_freq.upreg <-All_DEGs_Primary.Ambieant_v_Moderate %>%
            filter(Day_Trmt %in% c('Day7_Primary_AvM', 'Day14_Primary_AvM', 'Day21_Primary_AvM')) %>% 
            filter(up == TRUE) %>% 
            group_by(Row.names, up) %>%
            summarise(Freq = n()) %>% 
            arrange(desc(Freq))
DE_freq.upreg.sum <- DE_freq.upreg %>% group_by(Freq) %>% summarise(num.transcripts = n(),Direction = "upregulated")

DE_freq.downreg <-All_DEGs_Primary.Ambieant_v_Moderate %>%
            filter(Day_Trmt %in% c('Day7_Primary_AvM', 'Day14_Primary_AvM', 'Day21_Primary_AvM')) %>% 
            filter(down == TRUE) %>% 
            group_by(Row.names, down) %>%
            summarise(Freq = n()) %>%
            arrange(desc(Freq))
DE_freq.downreg.sum <- DE_freq.downreg %>% group_by(Freq) %>% summarise(num.transcripts = n(),Direction = "downregulated")

Constant.upreg.IDs <- DE_freq.upreg  %>%  dplyr::filter(Freq == 3)
Constant.downreg.IDs <- DE_freq.downreg  %>%  dplyr::filter(Freq == 3)

knitr::kable(Primary.DEGs.table, caption = "Primary Treatment DEGs [10 cpm in 50%; FDR 0.05; padj <0.05; LFC >1]")
knitr::kable(rbind(DE_freq.downreg.sum, DE_freq.upreg.sum),caption = "Primary Treatment DEGs, occurances on days 7, 14, 21: (1 = unique; 2 = d days; 3 = all days)")

write.csv((rbind(DE_freq.downreg, DE_freq.upreg)), "../Output/DESeq2/10cpm/Summary_Primary_DEGs_occurances_10CPM.csv")
write.csv((rbind(Constant.downreg.IDs, Constant.upreg.IDs)), "../Output/DESeq2/10cpm/Summary_Primary_DEGs_occurances_10CPM_GeneIDs.csv")
write.csv(Primary.DEGs.table, "../Output/DESeq2/10cpm/Summary_Primary_DEGs_10CPM.csv")
```


# Interaction terms: 
For loops using ~group -1 dds design to constrast all pairwise interaction terms: padj < 0.05; FDR 5%; LFC >1 (<-1) 
```{r GROUP_Day 7, 14, and 21}
# Day 7: ALL PAIRWISE COMPARISONS
d7.GROUP.DEGs_table = data.frame() # start a data frame
df_d7_group<- data.frame(resultsNames(dds.d7.group))
unique.vars <- df_d7_group[1,]
for(i in 1:nrow(df_d7_group)) {
  
      var <- df_d7_group[i,]
      T_F <- (df_d7_group[,1]!=c(unique.vars,var))
      df_d7_group_2 <- data.frame(df_d7_group[T_F,])

      for(j in 1:nrow(df_d7_group_2)) {
          
          var2 <- df_d7_group_2[j,]
          res <- results(dds.d7.group, contrast = list((var),(var2)),  alpha= 0.05) 
          res.table<-as.data.frame(table(res$padj<0.05))
          DEGs_total <- res.table[2,2]
          DEGs_total[is.na(DEGs_total)] <- 0
          DEGs.table <- data.frame(matrix(nrow = 1, ncol = 6)) # create a new data table
          colnames(DEGs.table)<-c('Day', 'Treatment', 'Pairwise.compare', 'DEGs_total', 'DEGs_UP', 'DEGs_DOWN') 
          res.order <- res[order(res$padj), ] ## Order by adjusted p-value
          DEGs.table$Day <- "7"
          DEGs.table$Treatment <- "full_interaction"
          DEGs.table$Pairwise.compare <- paste((substr(var, 14,15)), (substr(var2, 14,15)), sep = "_")
          DEGs.table$DEGs_total <- DEGs_total
          DEGs.table$DEGs_UP <- sum((res.order$log2FoldChange[1:(DEGs_total)] > 0) == TRUE) 
          DEGs.table$DEGs_DOWN <- sum((res.order$log2FoldChange[1:(DEGs_total)] < 0) == TRUE) 
          df.group <- data.frame(DEGs.table) # name dataframe for this single row
          d7.GROUP.DEGs_table <- rbind(d7.GROUP.DEGs_table, df.group) # bind to a cumulative list dataframe
          unique.vars <- unique((substr(d7.GROUP.DEGs_table$Pairwise.compare,1,2)))
      } 
}
d7.full.interaction.DEGs <- d7.GROUP.DEGs_table[-grep("MM_A|MA_A|MS_A|MM_M|MS_MA|AM_AA|AS_AA", d7.GROUP.DEGs_table$Pairwise.compare),] # ommit redundant pairwise results - same pairwise in opposite direaction

# Day 14: ALL PAIRWISE COMPARISONS
d14.GROUP.DEGs_table = data.frame() # start a data frame
df_d14_group<- data.frame(resultsNames(dds.d14.group))
unique.vars <- df_d14_group[1,]
for(i in 1:nrow(df_d14_group)) {
  
      var <- df_d14_group[i,]
      T_F <- (df_d14_group[,1]!=c(unique.vars,var))
      df_d14_group_2 <- data.frame(df_d14_group[T_F,])

      for(j in 1:nrow(df_d14_group_2)) {
          
          var2 <- df_d14_group_2[j,]
          res <- results(dds.d14.group, contrast = list((var),(var2)),  alpha= 0.05) 
          res.table<-as.data.frame(table(res$padj<0.05))
          DEGs_total <- res.table[2,2]
          DEGs_total[is.na(DEGs_total)] <- 0
          DEGs.table <- data.frame(matrix(nrow = 1, ncol = 6)) # create a new data table
          colnames(DEGs.table)<-c('Day', 'Treatment', 'Pairwise.compare', 'DEGs_total', 'DEGs_UP', 'DEGs_DOWN') 
          res.order <- res[order(res$padj), ] ## Order by adjusted p-value
          DEGs.table$Day <- "14"
          DEGs.table$Treatment <- "full_interaction"
          DEGs.table$Pairwise.compare <- paste((substr(var, 14,15)), (substr(var2, 14,15)), sep = "_")
          DEGs.table$DEGs_total <- DEGs_total
          DEGs.table$DEGs_UP <- sum((res.order$log2FoldChange[1:(DEGs_total)] > 0) == TRUE) 
          DEGs.table$DEGs_DOWN <- sum((res.order$log2FoldChange[1:(DEGs_total)] < 0) == TRUE) 
          df.group <- data.frame(DEGs.table) # name dataframe for this single row
          d14.GROUP.DEGs_table <- rbind(d14.GROUP.DEGs_table, df.group) # bind to a cumulative list dataframe
          unique.vars <- unique((substr(d14.GROUP.DEGs_table$Pairwise.compare,1,2)))
      } 
}
d14.full.interaction.DEGs <- d14.GROUP.DEGs_table[-grep("MM_A|MA_A|MS_A|MM_M|MS_MA|AM_AA|AS_AA", d14.GROUP.DEGs_table$Pairwise.compare),] # ommit redundant pairwise results - same pairwise in opposite direaction


# Day 21: ALL PAIRWISE COMPARISONS
d21.GROUP.DEGs_table = data.frame() # start a data frame
df_d21_group<- data.frame(resultsNames(dds.d21.group))
unique.vars <- df_d21_group[1,]
for(i in 1:nrow(df_d21_group)) {
  
      var <- df_d21_group[i,]
      T_F <- (df_d21_group[,1]!=c(unique.vars,var))
      df_d21_group_2 <- data.frame(df_d21_group[T_F,])

      for(j in 1:nrow(df_d21_group_2)) {
          
          var2 <- df_d21_group_2[j,]
          res <- results(dds.d21.group, contrast = list((var),(var2)),  alpha= 0.05) 
          res.table<-as.data.frame(table(res$padj<0.05))
          DEGs_total <- res.table[2,2]
          DEGs_total[is.na(DEGs_total)] <- 0
          DEGs.table <- data.frame(matrix(nrow = 1, ncol = 6)) # create a new data table
          colnames(DEGs.table)<-c('Day', 'Treatment', 'Pairwise.compare', 'DEGs_total', 'DEGs_UP', 'DEGs_DOWN') 
          res.order <- res[order(res$padj), ] ## Order by adjusted p-value
          DEGs.table$Day <- "21"
          DEGs.table$Treatment <- "full_interaction"
          DEGs.table$Pairwise.compare <- paste((substr(var, 14,16)), (substr(var2, 14,16)), sep = "_")
          DEGs.table$DEGs_total <- DEGs_total
          DEGs.table$DEGs_UP <- sum((res.order$log2FoldChange[1:(DEGs_total)] > 0) == TRUE) 
          DEGs.table$DEGs_DOWN <- sum((res.order$log2FoldChange[1:(DEGs_total)] < 0) == TRUE) 
          df.group <- data.frame(DEGs.table) # name dataframe for this single row
          d21.GROUP.DEGs_table <- rbind(d21.GROUP.DEGs_table, df.group) # bind to a cumulative list dataframe
          unique.vars <- unique((substr(d21.GROUP.DEGs_table$Pairwise.compare,1,3)))
      } 
}
View(d21.GROUP.DEGs_table)
d21.full.interaction.DEGs <- d21.GROUP.DEGs_table[-grep("MAA_A|MAM_A|MMA_A|MMM_A|MSA_A|MSM_A|MMA_MA|MMM_M|MSA_MA|MSM_MA|MAM_MAA", d21.GROUP.DEGs_table$Pairwise.compare),] # ommit redundant pairwise results - same pairwise in opposite direaction
View(d21.full.interaction.DEGs)
# knitr table
knitr::kable(rbind(d7.full.interaction.DEGs, d14.full.interaction.DEGs, d21.full.interaction.DEGs), caption = "All pairwise DEGs in 'group' DESeq2 model")
Summary_group_interactions <- rbind(d7.full.interaction.DEGs, d14.full.interaction.DEGs, d21.full.interaction.DEGs)
# Write csv 
write.csv(Summary_group_interactions, "../Output/DESeq2/10cpm/Summary_group_interactions.csv")
```

