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

```{r library_calls, echo=T, results='hide', message=FALSE, warning=FALSE}
library(tidyverse)
library(tximport) # Bioconductor
library(rhdf5) # Bioconductor
library(DESeq2) # Bioconductor
library(pheatmap)
library(ggthemes)
theme_set(theme_bw(base_size = 12))
```

```{r tximport, echo=T, results='hide', message=FALSE, warning=FALSE}
kallisto_dir <- "kallisto_results"
kallisto_df <- read_csv("kallisto_results/SraRunTable.csv")
row.names(kallisto_df) <- kallisto_df$Run

# Filter to include only one condition
coldata <- kallisto_df %>%
  filter(od600_nm == 2) %>%
  select(Run,
         od600_nm,
         phenotype)
# Put abundance file locations into a list
kallisto_files <- file.path(kallisto_dir,
                            coldata$Run,
                            "abundance.tsv")

kf2 <- file.path(kallisto_dir,
                            kallisto_df$Run,
                            "abundance.tsv")
txi <- tximport(kallisto_files, type="kallisto", txOut = TRUE)
txi2 <- tximport(kf2, type="kallisto", txOut = TRUE)

cts <- txi$counts # This is a count matrix like the one in the example
```

```{r DESeq, echo=T, results='hide', message=FALSE, warning=FALSE}
# Import tx object to DESeq format
dds <- DESeqDataSetFromTximport(txi,
                                colData = coldata,
                                design = ~ phenotype)

dds <- DESeq(dds) # This runs the standard DESeq normalization pipeline
resultsNames(dds) # Show contrast levels
res <- results(dds, contrast = c("phenotype", "HNS.KO", "WT")) # Get log2FC
results_df <- as.data.frame(res)
```

```{r counts, echo=T, results='hide', message=FALSE, warning=FALSE}
counts_df <- counts(dds) # Make a df with counts for each read
colnames(counts_df) <- coldata$Run # Fix column names
results_df <- cbind(counts_df, results_df)# Bind counts and results dfs
```

```{r plots, echo=T, results='show', message=FALSE, warning=FALSE}
select <- order(rowMeans(counts_df), decreasing=TRUE)[1:20] # Select top 20 genes
df <- as.data.frame(colData(dds)[,c("phenotype")]) #
rownames(df) <- coldata$Run # Fix rownames (must matach counts_df columns)
colnames(df) <- "phenotype" # Fix column name
# Make a heatmap
pheatmap(counts_df[select,],
         cluster_rows=FALSE,
         show_rownames=FALSE,
         cluster_cols=FALSE, annotation_col=df)

ntd <- normTransform(dds) # Normalizes to gives log2(n + 1)
colnames(ntd) <- coldata$Run # Fix column names
plotPCA(ntd, intgroup="phenotype") # Plot PCA
```