5_RNAseqVis_Scatter.Rmd

---
title: "Visualisation of DESeq2 results with scatter plots"
description: "DEG analysis based on DESeq2 and GSEA"
principal investigator: "Joaquín de Navascués"
researcher: "Joaquín de Navascués"
output: 
  html_document:
    toc: true
    toc_float: true
    code_folding: 'hide'
    theme: readable
    df_print: paged
---
```{r set-publication-theme, echo=FALSE, cache=FALSE}
ggplot2::theme_set(ggpubr::theme_pubr(base_size=10))
```

```{r setup, echo = FALSE, cache = FALSE}
knitr::opts_chunk$set(dev = 'png', 
                      fig.align = 'center', fig.height = 5, fig.width = 8.5, 
                      pdf.options(encoding = "ISOLatin9.enc"),
                      fig.path='scatter/figures/', warning=FALSE, message=FALSE)
```


**Libraries & external code:**
```{r load_libraries, warning=FALSE, echo=FALSE}
if (!require("librarian")) install.packages("librarian")
# data
librarian::shelf(
  # data manip.
  readxl, writexl, dplyr, tidyr, purrr, stringr, santoku, calibrate,
  # plotting
  cetcolor, ggplots2, ggpubr, ggthemes, ggtext, ggrepel,
# convenience
  here)

setwd(here::here()) # to distinguish from dplyr::here()

if(!exists("extract_regulated_sets2", mode="function")) source("utils.R")
```

**Path to definitive images (outside repo):**
```{r define_dir2figs}
figdir <- paste0(c(head(str_split(getwd(),'/')[[1]],-1),
                   paste0(tail(str_split(getwd(),'/')[[1]],1), '_figures')),
                 collapse='/')
dir.create(figdir, showWarnings = FALSE)
```


# RNAseq results visualisation with scatter plots


## 1 Getting ready


### 1.1 Load the DGE data

This gets us the DGE data from `DESeq2`, identified by FlyBase/Ensembl ID and gene symbol:
```{r load_DEG_data}
# experimental design and labels
targets <- readRDS('output/targets.RDS')
# DEG data
DaKD_deg <- readRDS('output/Control_vs_DaKD.RDS')
DaOE_deg <- readRDS('output/Control_vs_DaOE.RDS')
DaDaOE_deg <- readRDS('output/Control_vs_DaDaOE.RDS')
ScOE_deg <- readRDS('output/Control_vs_ScOE.RDS')
head(ScOE_deg, 1)
```

As discussed previously, we do not really need TPM values here, as the `baseMean` from the `DESeq2` normalisation does the work just as well and is more 'internally consistent'.

Now make a data frame with the status of each gene (up/down/non-regulated) for each condition:
```{r extract_regulated_sets2, warning=FALSE}
# classify genes as up-, down- or non-regulated by padj and fold change threshold
list_of_degs <- list('kd'=DaKD_deg, 'da'=DaOE_deg,
                     'dada'=DaDaOE_deg, 'sc'=ScOE_deg)
names_degs <- list('DaKD','DaOE','DaDaOE','ScOE')
cols.oi <- c('log2FoldChange', 'padj', 'baseMean',
             'gene_symbol', 'ensemblGeneID')
fc_thresh <- 1.5
regs <- extract_regulated_sets2(list_of_degs, names_degs, fc_thresh=fc_thresh, cols=cols.oi)

# attach the results to the DEG data frame for each condition
DaKD_deg$reg <- regs$DaKD
DaOE_deg$reg <- regs$DaOE
DaDaOE_deg$reg <- regs$DaDaOE
ScOE_deg$reg <- regs$ScOE
```

### 1.2 Select genes of interest

To create a list of genes of interest to track their individual expression in the different conditions, I am going to use a list of `genes_of_interest`, manually curated by JdN from the literature (including the Fly Cell Atlas paper) for a simplified cell type classification (ISC, EB, EE, EC, focused on the posterior midgut) and the markers defined by the Fly Cell Atlas consortium. Then I will remove redundancy and simplify the cell type classification.

- First make our list.

```{r load_GOI, message=FALSE}
genes_of_interest <- read_excel('resources/genes_of_interest.xlsx', sheet='core') %>%
  dplyr::select(gene_symbol, bonafide_celltype) %>%
  filter(bonafide_celltype != 'NA') %>%
  rename(celltype = bonafide_celltype) %>%
  group_by(celltype) %>%
  summarise(gene_symbol=paste(gene_symbol, collapse=', '))
```

- Then the FCA markers.

```{r loadFCA, message=FALSE}
markers <- read_xlsx(file.path(getwd(), 'resources', 'science.abk2432_table s2.xlsx')) %>%
  # get only the gut epithelium/muscle cell types
  filter(Tissue=='gut' | `Annotation...1`=='enteroendocrine cell') %>%
  # remove a hybrid/transitional cell type
  filter(!str_detect(`Annotation...1`, 'differentiating')) %>%
  # remove muscle/artefact cells
  filter(!str_detect(`Broad annotation`, 'muscle|artefact')) %>%
  # get just the columns with marker names
  dplyr::select(contains(c('Annotation...1', 'Marker'))) %>%
  # join all markers in a new column
  unite(gene_symbol, contains('Markers'), sep = ', ')
names(markers)[[1]] <- 'celltype'
```

- An merge the two:

```{r merge_markers, message=FALSE}
missd.spelld <- c(`wntD, `='', `LManIV, `='', `CG31269, `='', `CG43208, `='', `orcokinin B, `='',
                   AstB='Mip', poxn='Poxn')
markers <- rbind(genes_of_interest, markers) %>%
  # merge same cell type markers
  group_by(celltype) %>%
  summarise(gene_symbol=paste(gene_symbol, collapse=', ')) %>%
  # correct difference in synonym/capitalisation with DEG data...
  # ... while all gene symbols are in one string
  mutate(gene_symbol = str_replace_all(gene_symbol, missd.spelld)) %>%
  # make each gene name an independent string
  rowwise() %>% mutate(markers = list(unique(str_split(gene_symbol, ', ')[[1]])))
head(markers)
```

- Filter out the genes that do not have unique mapping to a cell type.

```{r filter_markers, message=FALSE}
# filter by exclusivity, but keeping what is common to all enterocytes
## first identify all EC markers
ECmarkers <- markers %>%
  filter(str_detect(celltype, 'enterocyte')) %>%
  dplyr::select(markers) %>%
  unlist(use.names = FALSE) %>%
  unique() %>% as.list()
## then non-EC cell markers, with all their repetitions!!
nonECmarkers <- markers %>%
  filter(!str_detect(celltype, 'enterocyte')) %>%
  dplyr::select(markers) %>%
  unlist(use.names = FALSE) %>% as.list()
## get all markers with one occurrence (≠unique!)
repmarkers <- unlist(c(ECmarkers, nonECmarkers))
reps <- rle(sort(repmarkers))
xmarkers <- reps$values[reps$lengths==1]
## apply criterion of exclusivity
exclude <- function(x) x[x %in% xmarkers]
markers <- markers %>%
    mutate(exclusive = list(exclude(markers)))
pmgcelltypes <- c('intestinal stem cell', 'enteroblast', 'enteroendocrine cell', 'enterocyte of posterior adult midgut epithelium')
pmg.markers <- markers %>%
  filter(celltype %in% pmgcelltypes) %>%
  dplyr::select(c(1,4)) %>%
  rename(xmarkers = exclusive)
```

- Add cell type colours HEX codes

```{r colour_markers}
write_xlsx(unnest_longer(pmg.markers, xmarkers) %>%
    rename(gene_symbol = xmarkers),
    path = 'output/preliminary_Table S4.xlsx')
cellcolours <- c('EB'='#56B2E9', 'EC'='#009E73', 'EE'='#CC79A7', 'ISC'='#D55E00')
pmg.markers$cellcolour <- cellcolours
```


## 2 MA plots using `ggpubr::ggmaplot` wrappers


MA plot for *daughterless* knockdown:
```{r ggmaplot2_da_knockdown}
repulsion <- list(box.padding = 0.1,
                  point.padding = 0.8,
                  nudge_x = 3,
                  nudge_y = 2,
                  force = 1,
                  force_pull = 0,
                  seed.up = 42,
                  seed.dn = 50,
                  xlims.up = c(14, NA),
                  ylims.up = c(5, NA),
                  xlims.dn = c(18, NA),
                  ylims.dn = c(NA, -1))

ggmaplot2(deg = DaKD_deg,
          fc_thresh = fc_thresh,
          markers = pmg.markers,
          md_label = '*esg > da^RNAi^*',
          repulsion = repulsion)
```

Now with cell type-specific colours:
```{r ggmaplot3_da_kd}
repulsion <- list(box.padding = 0.08,
                  point.padding = 0.1,
                  nudge_x = 0,
                  nudge_y = 0,
                  force = 1,
                  force_pull = -0.004,
                  seed.up = 42,
                  seed.dn = 5,
                  xlims.up = c(12, NA),
                  ylims.up = c(5.5, NA),
                  xlims.dn = c(16, NA),
                  ylims.dn = c(NA, -1))

ggmaplot3(deg = DaKD_deg,
          fc_thresh = fc_thresh,
          markers = pmg.markers, 
          md_label = '*esg > da^RNAi^*',
          repulsion = repulsion)


suppressMessages(ggsave('MA_daKD.pdf', plot = last_plot(), device = 'pdf',
                        path = figdir, dpi = 300))
```
  
```{r ggmaplot3_da_overexpression}
repulsion <- list(box.padding = 0.1,
                  point.padding = 0.1,
                  force = 1.5,
                  force_pull = -0.003,
                  nudge_x.up = 0,
                  nudge_y.up = 0,
                  nudge_x.dn = 0,
                  nudge_y.dn = 0,
                  seed.up = 42,
                  seed.dn = 5,
                  xlims.up = c(12, NA),
                  ylims.up = c(5, NA),
                  xlims.dn = c(16, NA),
                  ylims.dn = c(NA, -1))

ggmaplot3(deg = DaOE_deg,
          fc_thresh = fc_thresh,
          markers = pmg.markers, 
          md_label = '*esg > da*',
          repulsion = repulsion)

suppressMessages(ggsave('MA_daOE.pdf', plot = last_plot(), device = 'pdf',
                        path = figdir, dpi = 300))
```

To the Daughterless homodimer overexpression:
```{r ggmaplot2_da:da_overexpression}
repulsion <- list(box.padding = 0.01,
                  point.padding = 0.1,
                  nudge_x = 0,
                  nudge_y = 0,
                  force = 1,
                  force_pull = -0.004,
                  seed.up = 42,
                  seed.dn = 5,
                  xlims.up = c(16, NA),
                  ylims.up = c(3, NA),
                  xlims.dn = c(16, NA),
                  ylims.dn = c(NA, -2))

ggmaplot3(deg = DaDaOE_deg,
          fc_thresh = fc_thresh,
          markers = pmg.markers, 
          md_label = '*esg > da:da*',
          repulsion = repulsion)

suppressMessages(ggsave('MA_dadaOE.pdf', plot = last_plot(), device = 'pdf',
                        path = figdir, dpi = 300))
```

And, finally, to the *scute* overexpression:
```{r ggmaplot2_sc_overexpression}
repulsion <- list(box.padding = 0.01,
                  point.padding = 0.1,
                  nudge_x = 0,
                  nudge_y = 0,
                  force = 0.7,
                  force_pull = -0.003,
                  seed.up = 42,
                  seed.dn = 5,
                  xlims.up = c(8, NA),
                  ylims.up = c(7, NA),
                  xlims.dn = c(8, NA),
                  ylims.dn = c(NA, -6))

ggmaplot3(deg = ScOE_deg,
          fc_thresh = 2,
          markers = pmg.markers,
          md_label = '*esg > scute*',
          repulsion = repulsion)

suppressMessages(ggsave('MA_ScOE.pdf', plot = last_plot(), device = 'pdf',
                        path = figdir, dpi = 300))
```


## 3 Heatmap of individual genes


```{r}
tidy.markers <- pmg.markers %>% unnest_longer(xmarkers) %>%
  rename(gene_symbol = xmarkers)

DaKD_deg$condition <- '*da^RNAi^*'
DaOE_deg$condition <- '*da*'
DaDaOE_deg$condition <- '*da:da*'
ScOE_deg$condition <- '*scute*'
genehm.df <- bind_rows(ScOE_deg, DaDaOE_deg, DaOE_deg, DaKD_deg) %>%
  dplyr::select(ensemblGeneID, gene_symbol, log2FoldChange, padj, condition) %>%
  rename(p.adjust = padj) %>%
  filter(gene_symbol %in% tidy.markers$gene_symbol) %>%
  merge(., tidy.markers, by = 'gene_symbol') %>%
  mutate(condition = factor(condition, levels=c("*scute*", "*da:da*", "*da*", "*da^RNAi^*")))
```


```{r, fig.width=20}
p <- layer.heatmaph(genehm.df, arr='celltype')
p + geom_hline(aes(yintercept=3.5), linewidth = 0.5) +
  theme(plot.margin = margin(r = 25))
```

```{r, fig.width=20}
p <- layer.heatmaph(genehm.df, cluster=TRUE)
p + geom_hline(aes(yintercept=3.5), linewidth = 0.5) +
  theme(plot.margin = margin(r = 25))
```


```{r, fig.height = 20}
p <- layer.heatmapv(genehm.df, cluster=TRUE)
p + geom_vline(aes(xintercept=3.5), linewidth = 0.5)
suppressMessages(ggsave('marker_heatmap.pdf', plot = last_plot(), device = 'pdf',
                        path = figdir, dpi = 300))
```

```{r}
Espl.df <- bind_rows(ScOE_deg, DaDaOE_deg, DaOE_deg, DaKD_deg) %>%
  dplyr::select(ensemblGeneID, gene_symbol, log2FoldChange, padj, condition) %>%
  rename(p.adjust = padj) %>%
  filter(grepl("^E\\(spl\\)", gene_symbol)) %>%
  mutate(condition = factor(condition, levels=c("*scute*", "*da:da*", "*da*", "*da^RNAi^*"))) %>%
  mutate(cellcolour = "#000000")

p <- layer.heatmapv(Espl.df, cluster=TRUE)
p + geom_vline(aes(xintercept=3.5), linewidth = 0.5)
suppressMessages(ggsave('Esplit_heatmap.pdf', plot = last_plot(), device = 'pdf',
                        path = figdir, dpi = 300))
```