---
title: "Snow Crab Differential Gene Expression"
author: "Courtney Skalley"
format: revealjs
editor: visual
incremental: true
---

```{r, echo=FALSE, message=FALSE}
library(DESeq2)
library(tidyverse)
library(pheatmap)
library(RColorBrewer)
library(data.table)
```

```{r, cache=TRUE, echo=FALSE, message=FALSE}
countmatrix <- read.delim("../output/kallisto_01.isoform.counts.matrix", header = TRUE, sep = '\t')
rownames(countmatrix) <- countmatrix$X
countmatrix <- countmatrix[,-1]
```

```{r, cache=TRUE, echo=FALSE, message=FALSE}
countmatrix <- round(countmatrix, 0)
```

```{r, cache=TRUE, echo=FALSE, message=FALSE}
deseq2.colData <- data.frame(condition=factor(c(rep("control", 8), rep("low pH", 6))), 
                             type=factor(rep("single-read", 14)))
rownames(deseq2.colData) <- colnames(data)
deseq2.dds <- DESeqDataSetFromMatrix(countData = countmatrix,
                                     colData = deseq2.colData, 
                                     design = ~ condition)
dim(countmatrix)
dim(deseq2.colData)
deseq2.colData
deseq2.dds
```

```{r, cache=TRUE, echo=FALSE, message=FALSE}
deseq2.dds <- DESeq(deseq2.dds)
deseq2.res <- results(deseq2.dds)
deseq2.res <- deseq2.res[order(rownames(deseq2.res)), ]
```

## 

![Snow Crab (*Chionoecetes opilio*)](images/opilio_photo.jpeg){width="700"}

## Goal: Understand Snow Crab Response to Ocean Acidification

-   Compare snow crab response to ocean acidification
-   Conduct a differential gene expression analysis
-   Using data from Laura Spencer

## Experiment {.smaller}

63 juvenile snow crabs were exposed to different pH treatments for either a short (8 hours) or long (12 weeks) duration:

-   control (pH 8.1)
-   pH 7.8
-   pH 7.5

## RNA Seq Data {.smaller}

-   mRNA sequencing
-   Each of the 63 samples were run in both lanes
-   126 paired-end RNA-Seq data sets

## Available data

```{bash}
#list file names
ls /home/shared/8TB_HDD_01/snow_crab/5010/exp01
```

## Available data

```{bash}
#look at file
head -10 /home/shared/8TB_HDD_01/snow_crab/5010/exp01/5010_1_S1_L003_R1_001.fastq.gz
```


## Abundance estimates {.smaller}

-   Generate abundance estimates using Kallisto
-   Used  Trinity RNA-seq to assemble matrix to generate gene expression matrix  


## Abundance data {.smaller}

|    Contig    | X5010_1_S1_L003_R1 | X5010_1_S1_L003_R2 | X5010_2_S2_L003_R1 |
|--------------|:------------------:|:------------------:|:------------------:|
| contig_2632  | 16.2705            | 25.966000          |  6.88181           |
| contig_4568  | 0.0000             | 0.000000           |  0.00000           |
| contig_2580  | 300.0500           | 268.228000         |  212.03100         |
| contig_1896  | 2.0000             | 2.000000           |  0.00000           |
| contig_590   | 1007.6700          | 543.944000         |  506.597001        |
| contig_3193  | 25.7476            | 0.397017           |  7.17347           |


: Abundance matrix (estimates generated with Kallisto, matrix generated by Trinity RNA-seq assembly) 

## Code for Differential Gene Expression

``` {.r code-line-numbers="1-2"}
deseq2.colData <- data.frame(condition=factor(c(rep("control", 8), rep("low pH", 6))), 
                             type=factor(rep("single-read", 14)))
rownames(deseq2.colData) <- colnames(data)
deseq2.dds <- DESeqDataSetFromMatrix(countData = countmatrix,
                                     colData = deseq2.colData, 
                                     design = ~ condition)
dim(countmatrix)
dim(deseq2.colData)
deseq2.colData
deseq2.dds
```

## Code for Differential Gene Expression

``` {.r code-line-numbers="4-6"}
deseq2.colData <- data.frame(condition=factor(c(rep("control", 8), rep("low pH", 6))), 
                             type=factor(rep("single-read", 14)))
rownames(deseq2.colData) <- colnames(data)
deseq2.dds <- DESeqDataSetFromMatrix(countData = countmatrix,
                                     colData = deseq2.colData, 
                                     design = ~ condition)
dim(countmatrix)
dim(deseq2.colData)
deseq2.colData
deseq2.dds
```

## The Result: DGE for Snow Crabs

```{r}
#| echo: false
#| output-location: slide
tmp <- deseq2.res
# The main plot
plot(tmp$baseMean, tmp$log2FoldChange, pch=20, cex=0.45, ylim=c(-3, 3), log="x", col="darkgray",
     main="DEG pH 7.8  (pval <= 0.05)",
     xlab="mean of normalized counts",
     ylab="Log2 Fold Change")
# Getting the significant points and plotting them again so they're a different color
tmp.sig <- deseq2.res[!is.na(deseq2.res$padj) & deseq2.res$padj <= 0.05, ]
points(tmp.sig$baseMean, tmp.sig$log2FoldChange, pch=20, cex=0.45, col="red")
# 2 FC lines
abline(h=c(-1,1), col="blue")
```

## Next Steps

-   Conduct DGE using all 63 samples
-   Update code to combine forward and reverse sequences in DESeq
-   Visualize data in other ways
    -   Hierarchical Clustering Heatmap
    -   PCA
-   Create and annotate a de novo snow crab transcriptome
