differential-expression.qmd

---
title: "smallRNA Differential Expression"
author: "Lance Parsons"
date: last-modified
format:
  html:
    toc: true
    code-fold: true
    df-print: paged
    embed-resources: true
    page-layout: full
    fig-format: png
editor: source
---


# Load libraries

This project uses [`renv`](https://rstudio.github.io/renv/articles/renv.html)
to keep track of installed packages. Install `renv` if not installed and load
dependencies with `renv::restore()`.

```r
install.packages("renv")
renv::restore()
```

```{r}
#| label: load-packages
#| include: false
#| message: false
library("readr")
library("dplyr")
library("tidyr")
library("stringr")
library("Rsamtools")
library("GenomicAlignments")
library("plyranges")
library("tibble")
library("pheatmap")
library("ggplot2")
library("ggrepel")
library("ggnewscale")
library("svglite")
library("DESeq2")
library("apeglm")
library("RColorBrewer")
library("BiocParallel")
register(MulticoreParam(4))
```

# Read sample data

{{< include _sample-metadata.qmd >}}

# Exogenous RNA counts

{{< include _exogenous-alignments.qmd >}}

## Export exogenous RNA counts
```{r}
#| label: export-exogenous-rna-counts
# Split sample name into parts
exogenous_counts_temp <- exogenous_rna_count_data %>%
  mutate(sample_unit = str_replace_all(
    sample_unit,
    "_sorted|_pegRNA|_epegRNA",
    ""
  )) %>%
  tidyr::separate(sample_unit, "_", into = c(
    "background",
    "mixture",
    "day",
    "replicate"
  ), remove = FALSE)

# Calculate totals
exogenous_gene_count_totals <- exogenous_counts_temp %>%
  group_by(background, mixture, day, replicate, sequence_name, sample_unit) %>%
  summarize(count = sum(count), .groups = "keep") %>%
  mutate(category = "total", .before = count)

exogenous_counts_export <- rbind(
  exogenous_counts_temp,
  exogenous_gene_count_totals
) %>%
  dplyr::arrange(background, mixture, day, replicate, sequence_name, category)

# Export
write_tsv(exogenous_counts_export,
  file = "exexogenous_rna_count_data.tsv"
)
exogenous_counts_export
```

# Human small RNA gene counts

{{< include _human-small-rna-counts.qmd >}}

## Combine human and exogenous counts
```{r}
#| label: combine-human-exogenous-counts

# Reformat exogenous counts and filter out "other"
exogenous_gene_counts_category <- exogenous_rna_count_data %>%
  filter(category != "inactive_other") %>%
  tidyr::unite(gene, c(sequence_name, category), sep = ":") %>%
  tidyr::spread(sample_unit, count)

# Calculate and reformat exogenous total counts
exogenous_gene_counts_total <- exogenous_rna_mapped_totals %>%
  dplyr::ungroup() %>%
  dplyr::mutate(gene = paste(sequence_name, "total", sep = ":")) %>%
  dplyr::select(-"sequence_name") %>%
  tidyr::spread(sample_unit, mapped_fragments)

# Combine human and exogenous counts
gene_counts <- rbind(
  human_gene_counts,
  exogenous_gene_counts_category
)

gene_counts
```

## Exogenous RNA "gene" names
```{r}
#| label: list-exogenous-genes
# List of exogenous genes to highlight
exogenous_rna_names <- gene_counts %>%
  dplyr::filter(str_detect(gene, "^PJY")) %>%
  pull(gene)
```

## Import sample metadata and counts into DESeq2

Read these counts into DESeq2 along with the sample metadata.

Set the design to `~ exogenous_rna + day + cell_line`.

```{r}
#| label: create-deseq-dataset
dds <- DESeqDataSetFromMatrix(
  countData = gene_counts %>%
    tibble::column_to_rownames("gene") %>%
    replace(is.na(.), 0),
  colData = sample_units,
  design = ~ exogenous_rna + day + cell_line
)
dds
```

## Sample comparisons

### Variance Stabalized Transformation

```{r}
#| label: calculate-vst
vsd <- vst(dds, blind = FALSE)
```

### Heatmap of sample-to-sample distances

```{r}
#| label: sample-distance-matrix
sample_dists <- dist(t(assay(vsd)))

sample_dist_matrix <- as.matrix(sample_dists)
rownames(sample_dist_matrix) <- paste(vsd$cell_line,
  vsd$exogenous_rna,
  vsd$day,
  sep = "-"
)
colnames(sample_dist_matrix) <- NULL
```

```{r fig.height=8, fig.width=10}
#| label: sample-to-sample-heatmap
#| fig-height: 8
#| fig-width: 10
colors <- colorRampPalette(rev(brewer.pal(9, "Blues")))(255)
pheatmap(sample_dist_matrix,
  clustering_distance_rows = sample_dists,
  clustering_distance_cols = sample_dists,
  col = colors
)
```

### PCA plot of samples
```{r fig.height=4, fig.width=8}
#| label: sample-to-sample-pca-plot
#| fig-height: 4
#| fig-width: 8
plotPCA(vsd, intgroup = c("cell_line", "exogenous_rna", "day"))
```

# Differential Expression

## Plot aesthetics

{{< include _plot-aesthetics.qmd >}}

## Batch 2 - Day 1

Select subset of samples then calculate the log2 fold change:
`cell line` `P1E10` vs `Parental` 
```{r}
#| label: batch-2-day-1-diffexp
#| cache: true
dds_batch2_day1 <- subset(dds, select = colData(dds)$batch == "batch2" &
  colData(dds)$day == "day1")
dds_batch2_day1$exogenous_rna <- droplevels(dds_batch2_day1$exogenous_rna)
dds_batch2_day1$day <- droplevels(dds_batch2_day1$day)
design(dds_batch2_day1) <- ~ exogenous_rna + cell_line
dds_batch2_day1 <- DESeq(dds_batch2_day1, parallel = TRUE)
dds_batch2_day1
```

```{r}
#| label: batch-2-day-1-results
res_batch2_day1 <- results(
  dds_batch2_day1,
  name = "cell_line_P1E10_vs_Parental",
  alpha = 0.05
)
summary(res_batch2_day1)
```

### Log fold change shrinkage

Run [log fold change shrinkage](http://bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#log-fold-change-shrinkage-for-visualization-and-ranking) to account for low counts.

```{r}
#| label: batch-2-day-1-lfc-shrinkage
res_batch2_day1_lfc <- lfcShrink(
  dds_batch2_day1,
  coef = "cell_line_P1E10_vs_Parental",
  type = "apeglm",
  res = res_batch2_day1,
  parallel = TRUE
)
summary(res_batch2_day1_lfc)
```

### MA Plot (lfc shrunk)

MA Plot of shrunken log2 fold changes
```{r fig.height=8, fig.width=10}
#| label: batch-2-day-1-ma-plot
#| fig-height: 8
#| fig-width: 10
res_batch2_day1_lfc_df <- res_batch2_day1_lfc %>%
  as.data.frame() %>%
  rownames_to_column("gene") %>%
  dplyr::mutate(significant = as.factor(padj < 0.05)) %>%
  tidyr::separate_wider_delim(
    gene, ":",
    names = c("tmpgene", "category"),
    cols_remove = FALSE, too_few = "align_start"
  ) %>%
  left_join(
    rna_mixes %>% select(rna_species, rna_type),
    by = c("tmpgene" = "rna_species")
  ) %>%
  select(-tmpgene, -category)

write_tsv(res_batch2_day1_lfc_df, file = "diffexp_results_batch2_day1.tsv")

res_batch2_day1_lfc_labelled <- res_batch2_day1_lfc_df %>%
  dplyr::filter(gene %in% exogenous_rna_names) %>%
  na.omit() %>%
  separate(gene,
    into = c("rna_species", "category"),
    sep = ":",
    remove = FALSE
  ) %>%
  mutate(category = factor(category, levels = category_level_order)) %>%
  arrange(gene)

p <- ggplot(
  res_batch2_day1_lfc_df,
  aes(x = log2(baseMean), y = log2FoldChange, colour = significant)
) +
  scale_color_manual(
    name = "p-value",
    labels = sig_labels,
    values = human_gene_colors,
    na.value = na_color,
    drop = FALSE
  ) +
  scale_fill_manual(
    values = human_gene_colors,
    na.value = na_color,
    drop = FALSE
  ) +
  guides(colour = guide_legend(override.aes = list(size = 3))) +
  geom_point() +
  geom_hline(yintercept = 0, color = "black", lty = "dashed", lwd = 1) +
  ggrepel::geom_label_repel(
    data = res_batch2_day1_lfc_labelled,
    aes(
      label = rna_species,
      segment.colour = "black",
      fill = significant
    ),
    color = "white",
    min.segment.length = 0,
    max.overlaps = 30,
    size = 4,
    force_pull = 0.5,
    box.padding = 1,
    # Do not repel from top or bottom edges.
    ylim = c(-Inf, Inf),
    # Repel from 0 and from right edge
    xlim = c(0, NA),
    show.legend = FALSE
  ) +
  new_scale_color() +
  scale_color_manual(
    name = "Exogenous RNA Category",
    labels = category_labels,
    values = exogenous_gene_colors,
    na.value = na_color,
    drop = FALSE
  ) +
  geom_point(
    data = res_batch2_day1_lfc_labelled %>% arrange(category),
    aes(color = category),
    size = 4
  ) +
  xlab("log2 mean of normalized counts for all samples") +
  ylab("log2 fold change: cell line P1E10 vs Parental") +
  ggtitle("Batch 2 | Day 1") +
  theme_bw() +
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    axis.text = element_text(size = 24),
    axis.title = element_text(size = 24),
    plot.title = element_text(size = 24)
  )
print(p)
ggsave("diffexp-batch2-day1.svg", p)

# Remove labels
p_nolabels <- p
p_nolabels$layers <- p$layers[-3]
ggsave("diffexp-batch2-day1-nolabels.png", p_nolabels, dpi = 300)
```

### Counts of exogenous rna

Plot of counts for a single gene (with lowest adjusted p-value)
```{r fig.height=8, fig.width=15}
#| label: batch-2-day-1-exogenous-counts
#| fig-height: 8
#| fig-width: 15

plot_list <- list()
for (gene in res_batch2_day1_lfc_labelled$gene) {
  gene_prefix <- str_split(gene, "_")[[1]][1]
  d <- plotCounts(dds_batch2_day1,
    gene = gene,
    intgroup = "cell_line",
    returnData = TRUE
  )
  p <- ggplot(
    d %>%
      rownames_to_column("sample_unit") %>%
      dplyr::filter(grepl(gene_prefix, sample_unit)),
    aes(x = cell_line, y = count)
  ) +
    geom_point(position = position_jitter(width = 0.1, height = 0)) +
    ggtitle(gene)
  plot_list[[gene]] <- p
}
cowplot::plot_grid(plotlist = plot_list, ncol = 3)
```


## Batch 2 - Day 2

Select subset of samples then calculate the log2 fold change:
`cell line` `P1E10` vs `Parental` 
```{r}
#| label: batch-2-day-2-diffexp
#| cache: true
dds_batch2_day2 <- subset(dds, select = colData(dds)$batch == "batch2" &
  colData(dds)$day == "day2")
dds_batch2_day2$exogenous_rna <- droplevels(dds_batch2_day2$exogenous_rna)
dds_batch2_day2$day <- droplevels(dds_batch2_day2$day)
design(dds_batch2_day2) <- ~ exogenous_rna + cell_line
dds_batch2_day2 <- DESeq(dds_batch2_day2, parallel = TRUE)
dds_batch2_day2
```

```{r}
#| label: batch-2-day-2-results
res_batch2_day2 <- results(
  dds_batch2_day2,
  name = "cell_line_P1E10_vs_Parental",
  alpha = 0.05
)
summary(res_batch2_day2)
```

### Log fold change shrinkage

Run [log fold change shrinkage](http://bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#log-fold-change-shrinkage-for-visualization-and-ranking) to account for low counts.

```{r}
#| label: batch-2-day-2-lfc-shrinkage
res_batch2_day2_lfc <- lfcShrink(
  dds_batch2_day2,
  coef = "cell_line_P1E10_vs_Parental",
  type = "apeglm",
  res = res_batch2_day2,
  parallel = TRUE
)
summary(res_batch2_day2_lfc)
```

### MA Plot (lfc shrunk)

MA Plot of shrunken log2 fold changes
```{r fig.height=8, fig.width=10}
#| label: batch-2-day-2-ma-plot
#| fig-height: 8
#| fig-width: 10
res_batch2_day2_lfc_df <- res_batch2_day2_lfc %>%
  as.data.frame() %>%
  rownames_to_column("gene") %>%
  dplyr::mutate(significant = as.factor(padj < 0.05)) %>%
  tidyr::separate_wider_delim(
    gene, ":",
    names = c("tmpgene", "category"),
    cols_remove = FALSE, too_few = "align_start"
  ) %>%
  left_join(
    rna_mixes %>% select(rna_species, rna_type),
    by = c("tmpgene" = "rna_species")
  ) %>%
  select(-tmpgene, -category)

write_tsv(res_batch2_day2_lfc_df, file = "diffexp_results_batch2_day2.tsv")

res_batch2_day2_lfc_labelled <- res_batch2_day2_lfc_df %>%
  dplyr::filter(gene %in% exogenous_rna_names) %>%
  na.omit() %>%
  separate(gene,
    into = c("rna_species", "category"),
    sep = ":",
    remove = FALSE
  ) %>%
  mutate(category = factor(category, levels = category_level_order)) %>%
  arrange(gene)

p <- ggplot(
  res_batch2_day2_lfc_df,
  aes(x = log2(baseMean), y = log2FoldChange, color = significant)
) +
  scale_color_manual(
    name = "p-value",
    labels = sig_labels,
    values = human_gene_colors,
    na.value = na_color,
    drop = FALSE
  ) +
  scale_fill_manual(
    values = human_gene_colors,
    na.value = na_color,
    drop = FALSE
  ) +
  guides(colour = guide_legend(override.aes = list(size = 3))) +
  geom_point() +
  geom_hline(yintercept = 0, color = "black", lty = "dashed", lwd = 1) +
  ggrepel::geom_label_repel(
    data = res_batch2_day2_lfc_labelled,
    aes(
      label = rna_species,
      segment.colour = "black",
      fill = significant
    ),
    color = "white",
    min.segment.length = 0,
    max.overlaps = 30,
    size = 4,
    force_pull = 0.5,
    box.padding = 1,
    # Do not repel from top or bottom edges.
    ylim = c(-Inf, Inf),
    # Repel from 0 and from right edge
    xlim = c(0, NA),
    show.legend = FALSE
  ) +
  new_scale_color() +
  scale_color_manual(
    name = "Exogenous RNA Category",
    label = category_labels,
    values = exogenous_gene_colors,
    na.value = na_color,
    drop = FALSE
  ) +
  geom_point(
    data = res_batch2_day2_lfc_labelled %>% arrange(category),
    aes(color = category),
    size = 4
  ) +
  xlab("log2 mean of normalized counts for all samples") +
  ylab("log2 fold change: cell line P1E10 vs Parental") +
  ggtitle("Batch 2 | Day 2") +
  theme_bw() +
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    axis.text = element_text(size = 24),
    axis.title = element_text(size = 24),
    plot.title = element_text(size = 24)
  )
print(p)
ggsave("diffexp-batch2-day2.svg", p)

# Remove labels
p_nolabels <- p
p_nolabels$layers <- p$layers[-3]
ggsave("diffexp-batch2-day2-nolabels.png", p_nolabels, dpi = 300)
```

### Counts of exogenous rna

Plot of counts for a single gene (with lowest adjusted p-value)
```{r fig.height=8, fig.width=15}
#| label: batch-2-day-2-exogenous-counts
#| fig-height: 8
#| fig-width: 15

plot_list <- list()
for (gene in res_batch2_day2_lfc_labelled$gene) {
  gene_prefix <- str_split(gene, "_")[[1]][1]
  d <- plotCounts(dds_batch2_day2,
    gene = gene,
    intgroup = "cell_line",
    returnData = TRUE
  )
  p <- ggplot(
    d %>%
      rownames_to_column("sample_unit") %>%
      dplyr::filter(grepl(gene_prefix, sample_unit)),
    aes(x = cell_line, y = count)
  ) +
    geom_point(position = position_jitter(width = 0.1, height = 0)) +
    ggtitle(gene)
  plot_list[[gene]] <- p
}
cowplot::plot_grid(plotlist = plot_list, ncol = 3)
```


## Batch 3 - Day 1

Select subset of samples then calculate the log2 fold change:
`cell line` `P1E10` vs `Parental` 
```{r}
#| label: batch-3-day-1-diffexp
#| cache: true
dds_batch3_day1 <- subset(dds, select = colData(dds)$batch == "batch3" &
  colData(dds)$day == "day1")
dds_batch3_day1$exogenous_rna <- droplevels(dds_batch3_day1$exogenous_rna)
dds_batch3_day1$day <- droplevels(dds_batch3_day1$day)
design(dds_batch3_day1) <- ~ exogenous_rna + cell_line
dds_batch3_day1 <- DESeq(dds_batch3_day1, parallel = TRUE)
dds_batch3_day1
```

```{r}
#| label: batch-3-day-1-results
res_batch3_day1 <- results(
  dds_batch3_day1,
  name = "cell_line_P1E10_vs_Parental",
  alpha = 0.05
)
summary(res_batch3_day1)
```

### Log fold change shrinkage

Run [log fold change shrinkage](http://bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#log-fold-change-shrinkage-for-visualization-and-ranking) to account for low counts.


```{r}
#| label: batch-3-day-1-lfc-shrinkage
res_batch3_day1_lfc <- lfcShrink(
  dds_batch3_day1,
  coef = "cell_line_P1E10_vs_Parental",
  type = "apeglm",
  res = res_batch3_day1,
  parallel = TRUE
)
summary(res_batch3_day1_lfc)
```

### MA Plot (lfc shrunk)

MA Plot of shrunken log2 fold changes
```{r fig.height=8, fig.width=10}
#| label: batch-3-day-1-ma-plot
#| fig-height: 8
#| fig-width: 10
res_batch3_day1_lfc_df <- res_batch3_day1_lfc %>%
  as.data.frame() %>%
  rownames_to_column("gene") %>%
  dplyr::mutate(significant = as.factor(padj < 0.05)) %>%
  tidyr::separate_wider_delim(
    gene, ":",
    names = c("tmpgene", "category"),
    cols_remove = FALSE, too_few = "align_start"
  ) %>%
  left_join(
    rna_mixes %>% select(rna_species, rna_type) %>% unique(),
    by = c("tmpgene" = "rna_species")
  ) %>%
  select(-tmpgene, -category)

write_tsv(res_batch3_day1_lfc_df, file = "diffexp_results_batch3_day1.tsv")

res_batch3_day1_lfc_labelled <- res_batch3_day1_lfc_df %>%
  dplyr::filter(gene %in% exogenous_rna_names) %>%
  na.omit() %>%
  separate(gene,
    into = c("rna_species", "category"),
    sep = ":",
    remove = FALSE
  ) %>%
  mutate(category = factor(category, levels = category_level_order)) %>%
  arrange(gene)

p <- ggplot(
  res_batch3_day1_lfc_df,
  aes(x = log2(baseMean), y = log2FoldChange, colour = significant)
) +
  scale_color_manual(
    name = "p-value",
    labels = sig_labels,
    values = human_gene_colors,
    na.value = na_color,
    drop = FALSE
  ) +
  scale_fill_manual(
    values = human_gene_colors,
    na.value = na_color,
    drop = FALSE
  ) +
  guides(colour = guide_legend(override.aes = list(size = 3))) +
  geom_point() +
  geom_hline(yintercept = 0, color = "black", lty = "dashed", lwd = 1) +
  ggrepel::geom_label_repel(
    data = res_batch3_day1_lfc_labelled,
    aes(
      label = rna_species,
      segment.colour = "black",
      fill = significant
    ),
    color = "white",
    min.segment.length = 0,
    max.overlaps = Inf,
    size = 4,
    force_pull = 0.5,
    box.padding = 1,
    # Do not repel from top or bottom edges.
    ylim = c(-Inf, Inf),
    # Repel from 0 and from right edge
    xlim = c(0, NA),
    show.legend = FALSE
  ) +
  new_scale_color() +
  scale_color_manual(
    name = "Exogenous RNA Category",
    labels = category_labels,
    values = exogenous_gene_colors,
    na.value = na_color,
    drop = FALSE
  ) +
  geom_point(
    data = res_batch3_day1_lfc_labelled %>% arrange(category),
    aes(color = category),
    size = 4
  ) +
  xlab("log2 mean of normalized counts for all samples") +
  ylab("log2 fold change: cell line P1E10 vs Parental") +
  ggtitle("Batch 3 | Day 1") +
  theme_bw() +
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    axis.text = element_text(size = 24),
    axis.title = element_text(size = 24),
    plot.title = element_text(size = 24)
  )
print(p)
ggsave("diffexp-batch3-day1.svg", p)

# Remove labels
p_nolabels <- p
p_nolabels$layers <- p$layers[-3]
ggsave("diffexp-batch3-day1-nolabels.png", p_nolabels, dpi = 300)
```

### Counts of exogenous rna

Plot of counts for a single gene (with lowest adjusted p-value)
```{r fig.height=30, fig.width=15}
#| label: batch-3-day-1-exogenous-counts
#| fig-height: 30
#| fig-width: 15

res_batch3_day1_lfc_labelled <- res_batch3_day1_lfc_df %>%
  dplyr::filter(gene %in% exogenous_rna_names) %>%
  na.omit() %>%
  arrange(gene)

plot_list <- list()
for (gene in res_batch3_day1_lfc_labelled$gene) {
  gene_rna_type <- str_split(gene, ":")[[1]][1]
  exogenous_rna <- rna_mixes %>%
    dplyr::filter(rna_species == gene_rna_type) %>%
    pull(exogenous_rna)
  exogenous_rna_regex <- paste(exogenous_rna, collapse = "|")

  d <- plotCounts(dds_batch3_day1,
    gene = gene,
    intgroup = "cell_line",
    returnData = TRUE
  )
  p <- ggplot(
    d %>%
      rownames_to_column("sample_unit") %>%
      dplyr::filter(grepl(exogenous_rna_regex, sample_unit)),
    aes(x = cell_line, y = count)
  ) +
    geom_point(position = position_jitter(width = 0.1, height = 0)) +
    ggtitle(gene)
  plot_list[[gene]] <- p
}
cowplot::plot_grid(plotlist = plot_list, ncol = 3)
```



## Batch 3 - Day 2

Select subset of samples then calculate the log2 fold change:
`cell line` `P1E10` vs `Parental` 
```{r}
#| label: batch-3-day-2-diffexp
#| cache: true
dds_batch3_day2 <- subset(dds, select = colData(dds)$batch == "batch3" &
  colData(dds)$day == "day2")
dds_batch3_day2$exogenous_rna <- droplevels(dds_batch3_day2$exogenous_rna)
dds_batch3_day2$day <- droplevels(dds_batch3_day2$day)
design(dds_batch3_day2) <- ~ exogenous_rna + cell_line
dds_batch3_day2 <- DESeq(dds_batch3_day2, parallel = TRUE)
dds_batch3_day2
```

```{r}
#| label: batch-3-day-2-results
res_batch3_day2 <- results(
  dds_batch3_day2,
  name = "cell_line_P1E10_vs_Parental",
  alpha = 0.05
)
summary(res_batch3_day2)
```


### Log fold change shrinkage

Run [log fold change shrinkage](http://bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#log-fold-change-shrinkage-for-visualization-and-ranking) to account for low counts.


```{r}
#| label: batch-3-day-2-lfc-shrinkage
res_batch3_day2_lfc <- lfcShrink(
  dds_batch3_day2,
  coef = "cell_line_P1E10_vs_Parental",
  type = "apeglm",
  res = res_batch3_day2,
  parallel = TRUE
)
summary(res_batch3_day2_lfc)
```

### MA Plot (lfc shrunk)

MA Plot of shrunken log2 fold changes
```{r fig.height=8, fig.width=10}
#| label: batch-3-day-2-ma-plot
#| fig-height: 8
#| fig-width: 10
res_batch3_day2_lfc_df <- res_batch3_day2_lfc %>%
  as.data.frame() %>%
  rownames_to_column("gene") %>%
  dplyr::mutate(significant = as.factor(padj < 0.05)) %>%
  tidyr::separate_wider_delim(
    gene, ":",
    names = c("tmpgene", "category"),
    cols_remove = FALSE, too_few = "align_start"
  ) %>%
  left_join(
    rna_mixes %>% select(rna_species, rna_type) %>% unique(),
    by = c("tmpgene" = "rna_species")
  ) %>%
  select(-tmpgene, -category)

write_tsv(res_batch3_day2_lfc_df, file = "diffexp_results_batch3_day2.tsv")

res_batch3_day2_lfc_labelled <- res_batch3_day2_lfc_df %>%
  dplyr::filter(gene %in% exogenous_rna_names) %>%
  na.omit() %>%
  separate(gene,
    into = c("rna_species", "category"),
    sep = ":",
    remove = FALSE
  ) %>%
  mutate(category = factor(category, levels = category_level_order)) %>%
  arrange(gene)

p <- ggplot(
  res_batch3_day2_lfc_df,
  aes(x = log2(baseMean), y = log2FoldChange, colour = significant)
) +
  scale_color_manual(
    name = "p-value",
    labels = sig_labels,
    values = human_gene_colors,
    na.value = na_color,
    drop = FALSE
  ) +
  scale_fill_manual(
    values = human_gene_colors,
    na.value = na_color,
    drop = FALSE
  ) +
  guides(colour = guide_legend(override.aes = list(size = 3))) +
  geom_point() +
  geom_hline(yintercept = 0, color = "black", lty = "dashed", lwd = 1) +
  ggrepel::geom_label_repel(
    data = res_batch3_day2_lfc_labelled,
    aes(
      label = rna_species,
      segment.colour = "black",
      fill = significant
    ),
    color = "white",
    min.segment.length = 0,
    max.overlaps = Inf,
    size = 4,
    force_pull = 0.5,
    box.padding = 1,
    # Do not repel from top or bottom edges.
    ylim = c(-Inf, Inf),
    # Repel from 0 and from right edge
    xlim = c(0, NA),
    show.legend = FALSE
  ) +
  new_scale_color() +
  scale_color_manual(
    name = "Exogenous RNA Category",
    labels = category_labels,
    values = exogenous_gene_colors,
    na.value = na_color,
    drop = FALSE
  ) +
  geom_point(
    data = res_batch3_day2_lfc_labelled %>% arrange(category),
    aes(color = category),
    size = 4
  ) +
  xlab("log2 mean of normalized counts for all samples") +
  ylab("log2 fold change: cell line P1E10 vs Parental") +
  ggtitle("Batch 3 | Day 2") +
  theme_bw() +
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    axis.text = element_text(size = 24),
    axis.title = element_text(size = 24),
    plot.title = element_text(size = 24)
  )
print(p)
ggsave("diffexp-batch3-day2.svg", p)

# Remove labels
p_nolabels <- p
p_nolabels$layers <- p$layers[-3]
ggsave("diffexp-batch3-day2-nolabels.png", p_nolabels, dpi = 300)
```

### Counts of exogenous rna

Plot of counts for a single gene (with lowest adjusted p-value)
```{r fig.height=30, fig.width=15}
#| label: batch-3-day-2-exogenous-counts
#| fig-height: 30
#| fig-width: 15

res_batch3_day2_lfc_labelled <- res_batch3_day2_lfc_df %>%
  dplyr::filter(gene %in% exogenous_rna_names) %>%
  na.omit() %>%
  arrange(gene)

plot_list <- list()
for (gene in res_batch3_day2_lfc_labelled$gene) {
  gene_rna_type <- str_split(gene, ":")[[1]][1]
  exogenous_rna <- rna_mixes %>%
    dplyr::filter(rna_species == gene_rna_type) %>%
    pull(exogenous_rna)
  exogenous_rna_regex <- paste(exogenous_rna, collapse = "|")

  d <- plotCounts(dds_batch3_day2,
    gene = gene,
    intgroup = "cell_line",
    returnData = TRUE
  )
  p <- ggplot(
    d %>%
      rownames_to_column("sample_unit") %>%
      dplyr::filter(grepl(exogenous_rna_regex, sample_unit)),
    aes(x = cell_line, y = count)
  ) +
    geom_point(position = position_jitter(width = 0.1, height = 0)) +
    ggtitle(gene)
  plot_list[[gene]] <- p
}
cowplot::plot_grid(plotlist = plot_list, ncol = 3)
```