---
title: "Week 2 - Differential Gene Expression"
author: Haila Schultz
output: 
  html_document: 
    theme: readable
date: "2023-04-13"
---

This code uses RNA-Seq reads to examine differential gene expression for two treatments of an experiment.

# Load packages


Load R packages
```{r load packages, message=FALSE}
library(DESeq2)
library(tidyverse)
library(pheatmap)
library(RColorBrewer)
library(data.table)
library(DT)
```


# Kallisto alignment

make sure kallisto works
```{r check kallisto, eval=FALSE, engine="bash"}
#kallisto is already downloaded in Raven
/home/shared/kallisto/kallisto -h
```


download reference
```{r download reference, eval=FALSE, engine="bash"}
#move to preferred directory
cd ../data

#download from url
curl -O https://gannet.fish.washington.edu/seashell/bu-github/nb-2023/Cgigas/data/rna.fna
```


create index and rename
```{r create index, eval=FALSE, engine="bash"}
/home/shared/kallisto/kallisto \
index -i \
../data/cgigas_roslin_rna.index \
../data/rna.fna
```


Download data ending in ".fastq.gz". wget downloads multiple files recursively.
```{r download data, eval=FALSE, engine="bash"}
wget --recursive --no-parent --no-directories \
--no-check-certificate \
--accept '*.fastq.gz' \
https://gannet.fish.washington.edu/seashell/bu-github/nb-2023/Cgigas/data/nopp/
```


Create subdirectory and align reads with kallisto
```{r align reads, eval=FALSE, engine="bash"}
#create subdirectory
mkdir ../output/kallisto_01

#run kallisto
find ../data/*fastq.gz \
| xargs basename -s _L001_R1_001.fastq.gz | xargs -I{} /home/shared/kallisto/kallisto \
quant -i ../data/cgigas_roslin_rna.index \
-o ../output/kallisto_01/{} \
-t 40 \
--single -l 100 -s 10 ../data/{}_L001_R1_001.fastq.gz
```


Create gene expression matrix
```{r gene expression matrix, eval=FALSE, engine="bash"}
perl /home/shared/trinityrnaseq-v2.12.0/util/abundance_estimates_to_matrix.pl \
--est_method kallisto \
    --gene_trans_map none \
    --out_prefix ../output/kallisto_01 \
    --name_sample_by_basedir \
    ../output/kallisto_01/D54_S145/abundance.tsv \
    ../output/kallisto_01/D56_S136/abundance.tsv \
    ../output/kallisto_01/D58_S144/abundance.tsv \
    ../output/kallisto_01/M45_S140/abundance.tsv \
    ../output/kallisto_01/M48_S137/abundance.tsv \
    ../output/kallisto_01/M89_S138/abundance.tsv \
    ../output/kallisto_01/D55_S146/abundance.tsv \
    ../output/kallisto_01/D57_S143/abundance.tsv \
    ../output/kallisto_01/D59_S142/abundance.tsv \
    ../output/kallisto_01/M46_S141/abundance.tsv \
    ../output/kallisto_01/M49_S139/abundance.tsv \
    ../output/kallisto_01/M90_S147/abundance.tsv \
    ../output/kallisto_01/N48_S194/abundance.tsv \
    ../output/kallisto_01/N50_S187/abundance.tsv \
    ../output/kallisto_01/N52_S184/abundance.tsv \
    ../output/kallisto_01/N54_S193/abundance.tsv \
    ../output/kallisto_01/N56_S192/abundance.tsv \
    ../output/kallisto_01/N58_S195/abundance.tsv \
    ../output/kallisto_01/N49_S185/abundance.tsv \
    ../output/kallisto_01/N51_S186/abundance.tsv \
    ../output/kallisto_01/N53_S188/abundance.tsv \
    ../output/kallisto_01/N55_S190/abundance.tsv \
    ../output/kallisto_01/N57_S191/abundance.tsv \
    ../output/kallisto_01/N59_S189/abundance.tsv
```


# Run DESeq2

Now that the reads have been downloaded and aligned, we can conduct differential gene expression analysis.

Read in count matrix
```{r read in count matrix}
#read matrix into r
countmatrix <- read.delim("../output/kallisto_01.isoform.counts.matrix", header = TRUE, sep = '\t')
rownames(countmatrix) <- countmatrix$X
countmatrix <- countmatrix[,-1]

#explore table
head(countmatrix)
```


round integers up
```{r round integers}
countmatrix <- round(countmatrix, 0)
str(countmatrix)
```


Get DEGs based on Desication
```{r get DEGs, message=FALSE}
deseq2.colData <- data.frame(condition=factor(c(rep("control", 12), rep("desicated", 12))), 
                             type=factor(rep("single-read", 24)))
rownames(deseq2.colData) <- colnames(data)
deseq2.dds <- DESeqDataSetFromMatrix(countData = countmatrix,
                                     colData = deseq2.colData, 
                                     design = ~ condition)
```

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

Explore results 
```{r explore results}
head(deseq2.res)
```
# Visualize Data

Count number of hits with adjusted p-value less then 0.05
```{r count hits}
dim(deseq2.res[!is.na(deseq2.res$padj) & deseq2.res$padj <= 0.05, ])
```

plot significant points
```{r plot significant points, warning=FALSE}
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 Dessication  (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")
```

make table
```{r make table from desication}
write.table(tmp.sig, "../output/DEGlist.tab", sep = '\t', row.names = T)
deglist <- read.csv("../output/DEGlist.tab", sep = '\t', header = TRUE)
deglist$RowName <- rownames(deglist)
datatable(deglist)
```


Look at top 9 some genes of interest
```{r find top nine genes}
#examine data
colData(deseq2.dds)
# get results and arrange by pvalue
res <- results(deseq2.dds, tidy=TRUE, contrast=c("condition", "control", "desicated")) %>%
  arrange(padj, pvalue) %>%
  as_tibble()

# Define the top 10 genes of interest
goi <- res$row[1:9]
stopifnot(all(goi %in% names(deseq2.dds)))
goi
```
make count matrix
```{r top nine gens count matrix}
#make matrix
tcounts <- t(log2((counts(deseq2.dds[goi, ], normalized=TRUE, replaced=FALSE)+.5))) %>%
  merge(colData(deseq2.dds), ., by="row.names") %>%
  gather(gene, expression, (ncol(.)-length(goi)+1):ncol(.))

#take a look at matrix
tcounts %>% 
  select(Row.names, condition, gene, expression) %>% 
  head %>% 
  knitr::kable()
```

create box plot
```{r top nine genes boxplot}
ggplot(tcounts, aes(condition, expression)) + 
  geom_boxplot() + 
  facet_wrap(~gene, scales="free_y") + 
  labs(x="treatment", 
       y="Expression (log normalized counts)", 
       title="Top Nine Gene Results")+
  theme_bw()
```

