analyze_host_fold_changes.Rmd

---
title: "Determine list of induced genes under various treatments"
author: "`r Sys.info()[['user']]`"
date: "`r format(Sys.time(), '%d %B, %Y')`"
output:
  github_document:
    toc: yes
---

We look for genes that are induced across both Ribo-seq and RNA-seq samples.
Therefore, we treat these two measurements for each treatment as replicates.

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

# Import libraries and define analysis specific variables
```{r}
library(DESeq2)
library(tidyverse)

# set the lowest counts to be these for DESeq2 fold-changes if it is zero
deseq2_floor <- 1
```

# Read gene annotations

Select only genes that have annotated protein coding sequences in the CCDS 
database.
```{r}
gene_annotations <- read_tsv(glue::glue("ftp://ftp.ebi.ac.uk/pub/databases/",
  "genenames/new/tsv/locus_groups/protein-coding_gene.txt")) %>%
  filter(!is.na(ccds_id)) %>%
  select(ensembl_gene_id, symbol) %>%
  rename(gene_id = ensembl_gene_id) %>% 
  print()
```

# Read in gene counts and process it for input to DESeq2

Exclude Ribo-seq  + LTM data for fold-change calculation

```{r}
counts <- list.files("../processeddata/", pattern = "gencode.genes.results",
                          recursive = T,
                          full.names = T) %>% 
  enframe("sno", "file") %>% 
  # exclude LTM data
  filter(!str_detect(file, "ltm")) %>% 
  # get sample name
  mutate(sample = str_extract(file, "[^/]+(?=/gencode.genes.results)")) %>% 
  # read data in
  mutate(data = map(file, data.table::fread)) %>% 
  unnest(data) %>% 
  # remove unwanted columns
  select(-sno, -file, -length, -effective_length, -TPM, -FPKM, -`transcript_id(s)`) %>% 
  mutate(expected_count = as.integer(expected_count)) %>% 
  # first select only gene, sample pairs with non-zero counts
  filter(expected_count > 0) %>% 
  # next create all sample-gene pairs filling with deseq2_floor if it doesn't exist
  complete(sample, gene_id, fill = list("expected_count" = deseq2_floor)) %>% 
  # select only rows that have a minimum of 100 counts across all samples
  group_by(gene_id) %>% 
  mutate(total_counts = sum(expected_count)) %>% 
  ungroup() %>% 
  filter(total_counts > 100) %>% 
  select(-total_counts)

gene_counts <- counts %>% 
  # create one column per sample
  spread(sample, expected_count) %>% 
  # remove gene id suffix to match with gene_name
  mutate(gene_id = str_extract(gene_id, "\\w+"))

print(gene_counts)

gene_counts <- gene_counts %>% 
  as.data.frame() %>% 
  column_to_rownames("gene_id")
```

# Run DESeq2 to get log fold changes

Treat Ribo and mRNA as replicates to get non-zero P-values
```{r}
coldata <- colnames(gene_counts) %>% 
  enframe("sno", "sample") %>% 
  mutate(treatment = str_extract(sample, "(?<=mrna_|cyclo_).+$")) %>% 
  mutate(type = str_extract(sample, "^(mrna|cyclo)")) %>% 
  mutate(colname = sample) %>% 
  as.data.frame() %>% 
  column_to_rownames("colname") %>% 
  select(-sno) %>% 
  print()

ddsObject <- DESeqDataSetFromMatrix(countData = gene_counts,
                                    colData = coldata,
                                    design = ~ treatment)

dds <- DESeq(ddsObject)
```

# Get fold-changes between specific sample pairs

The sample pairs were created manually in `../tables/samplepairs_for_deseq2.tsv`.

```{r}
lfc <- read_tsv("../tables/samplepairs_for_deseq2.tsv") %>% 
  mutate(deseq_results = map2(treatment1, treatment2, function(x, y) 
    results(dds, contrast = c("treatment", x, y)))) %>%
  mutate(lfc = map(deseq_results, function(res)
    res %>% as.data.frame() %>% rownames_to_column("gene_id"))) %>%
  select(-deseq_results) %>%
  unnest() %>%
  print()
```

# Get list of inferferon-induced genes

We treat Ribo and mRNA as replicates
```{r}
lfc %>% 
  filter(treatment1 == "ifn") %>% 
  left_join(gene_annotations) %>% 
  rename(lfc = log2FoldChange) %>% 
  select(-stat, -lfcSE, -padj) %>% 
  filter(pvalue < 0.001 & lfc >= log2(2)) %>%
  arrange(pvalue) %>% 
  write_tsv("../tables/interferon_induced_genes.tsv") %>% 
  print()
```

# Get list of virus-induced genes 

We treat Ribo and mRNA as replicates

```{r}
lfc %>% 
  filter(treatment1 == "vir") %>% 
  left_join(gene_annotations) %>% 
  rename(lfc = log2FoldChange) %>% 
  select(-stat, -lfcSE, -padj) %>% 
  filter(pvalue < 0.001 & lfc >= log2(2)) %>%
  arrange(pvalue) %>% 
  write_tsv("../tables/virus_induced_genes.tsv") %>% 
  print()
```

# Get list of inferferon + virus induced genes

We treat Ribo and mRNA as replicates

```{r}
lfc %>% 
  filter(treatment1 == "ifn_vir") %>% 
  left_join(gene_annotations) %>% 
  rename(lfc = log2FoldChange) %>% 
  select(-stat, -lfcSE, -padj) %>% 
  filter(pvalue < 0.001 & lfc >= log2(2)) %>%
  arrange(pvalue) %>% 
  write_tsv("../tables/interferon_plus_virus_induced_genes.tsv") %>% 
  print()
```