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

# Load libraries

```{r}
library(DESeq2)
library(plyranges)
library(tidyverse)
```

# Load RNA-Seq gene counts 

```{r}
counts <- list.files('data/gene_counts/', pattern = '.tsv.gz$', full.names = T) %>% 
  as_tibble_col('file') %>% 
  mutate(data = map(file, read_tsv)) %>% 
  mutate(sample_name = str_extract(file, '[^/]+(?=.tsv.gz$)')) %>% 
  select(-file) %>%
  unnest() %>% 
  group_by(gene_name, sample_name) %>% 
  slice(1) %>% 
  ungroup() %>% 
  print()
```

# Prepare count matrix

```{r}
countdata <- counts %>% 
  pivot_wider(names_from = sample_name, values_from = count) %>%
  column_to_rownames("gene_name")

head(countdata)
```

# Prepare annotations matrix

```{r}
coldata <- colnames(countdata) %>% 
  as_tibble_col("sample_name") %>% 
  left_join(read_tsv("sample_annotations.tsv"), by = "sample_name") %>% 
  select(sample_name, shrna, replicate) %>%
  column_to_rownames("sample_name") %>%
  print()
```

# Run differential expression analysis

```{r}
ddsObject <- DESeq2::DESeqDataSetFromMatrix(countData = countdata,
                                            colData = coldata, design = ~ shrna)
dds <- DESeq2::DESeq(ddsObject)
```
# Extract differentially expressed genes and write to output

```{r}
lfc <- list(sample1 = c("rpl5", "rps19", "tsr2"), sample2 = c("luc", "luc", "luc")) %>% 
  as_tibble() %>%
  mutate(deseq_results = map2(sample1, sample2, function(x, y)
    DESeq2::results(dds, contrast = c("shrna", x, y)))) %>%
  mutate(lfc = map(deseq_results, function(res)
    res[c('baseMean', 'log2FoldChange', 'pvalue')] %>% as.data.frame() %>% rownames_to_column("gene_name"))) %>%
  select(-deseq_results) %>%
  unnest() %>%
  rename(mean_counts = baseMean, lfc = log2FoldChange) %>%
  filter(pvalue < 0.05) %>%
  arrange(lfc) %>%
  write_tsv("../data/differentially_expressed_genes.tsv") %>%
  print()
```
```

