---
title: "RNA sequencing data analysis"
author: "Cihan Erkut"
date: "20240109"
date-modified: "`r Sys.time()`"
format:
  html:
    theme: materia
    toc-title: Contents
    highlight-style: github
    df-print: tibble
    toc: true
    number-sections: true
    embed-resources: true
    code-line-numbers: true
    code-block-border-left: true
    code-fold: true
    code-tools: true
knitr:
  opts_chunk: 
    tidy: "styler"
execute:
  cache: true
---

## Setup

Load libraries and define custom functions.

```{r}
#| label: setup
#| message: false

# Load libraries in this order to avoid namespace collision.
# We also need ashr and IHW packages for DESeq.
# 
library(ComplexHeatmap) # BioConductor
library(DESeq2)         # BioConductor
library(sva)            # BioConductor
library(zFPKM)          # BioConductor
library(glue)           # CRAN
library(gplots)         # CRAN
library(tidyverse)      # CRAN
library(ggrepel)        # CRAN
library(GGally)         # CRAN
library(UpSetR)         # CRAN
library(ggthemes)       # CRAN
library(fs)             # CRAN
library(writexl)        # CRAN
library(here)           # CRAN

```

```{r}
#| label:  import functions

import::here("parse_featureCounts", .from = "functions.R")
import::here("select_expressed_genes", .from = "functions.R")
import::here("deseq_data", .from = "functions.R")
import::here("analyze_batch_effects", .from = "functions.R")
import::here("get_results", .from = "functions.R")

```

## Load data

```{r}
#| label: load data

# RNA sequencing library was stranded, therefore we only use reversely-stranded
# read count (num_reads_rv) and its derivatives (FPKM, TPM).
# 
# Sample information is extracted from path name and separated into treatment
# (e.g. A2, B1, D4, E35) and replicate (e.g. 4-6).
# 
featureCounts_files <- dir_ls(
  path = here("raw"),
  glob = "*.fpkm_tpm.featureCounts.tsv"
)

table_darpin <- map_dfr(
  .x = featureCounts_files,
  .f = parse_featureCounts,
  .id = "Path"
) |>
  separate_wider_regex(
    cols = Path,
    patterns = c(
      ".*/cell-",
      Treatment = ".*",
      "-",
      Replicate = "\\d?",
      "_B290.*"
    )
  ) |>
  unite(
    col = "Group",
    Treatment,
    Replicate,
    remove = FALSE
  ) |>
  mutate(across(Group:Replicate, as_factor)) |>
  rename(
    counts = num_reads_rv,
    FPKM = FPKM_rv,
    TPM = TPM_rv
  )

# Read GRCh37 Gencode annotations
table_annotation <- read_tsv(
  file = here("data/GRCh37_annotation.tsv")
)

```

## Data preparation

Expressed genes are selected using zFPKM method. A zFPKM cutoff of -3 is used as recommended by the developers of the method. Below are individual plots of Gaussian fits that are used to select expressed genes.

```{r}
#| label: zFPKM DARPin

# Choose reliably expressed genes for each dataset
#
genes_expressed_darpin <- select_expressed_genes(table_darpin)

```

In DARPin experiment there are **`r length(genes_expressed_darpin)`** genes that are reliably expressed.

```{r}
#| label: prepare for DESeq

# Prepare the count matrix and column data to be used in constructing the DESeq
# datasets
#
deseq_data_darpin <- deseq_data(
  table_expression = table_darpin,
  genes_expressed = genes_expressed_darpin
)

matrix_darpin_counts <- deseq_data_darpin$matrix_counts
table_darpin_colData <- deseq_data_darpin$table_colData

```

## Initial analysis

We search for batch effects in using PCA on rlog-transformed counts. It looks like in the DARPin experiment replicate 4 is separated from replicates 5 and 6 in all treatments except E35.

```{r}
#| label: batch darpin
#| message: false
#| fig-height: 4
#| fig-width: 5

# Batch effect analysis
# 
matrix_darpin_rlog_corr <- analyze_batch_effects(
  matrix_counts = matrix_darpin_counts,
  table_colData = table_darpin_colData,
  return_corrected_matrix = TRUE
)

```

Batch-corrected gene expression levels (i.e. rlog-transformed counts) of reliably-expressed genes were averaged for each gene in each treatment. Then we explore correlation between each treatment.

Correlation coefficients are significant and high, but the significance is probably misleading because of the very high number of data points.

```{r}
#| label: darpin correlation

# Correlation of gene expression among treatments
# 
matrix_darpin_rlog_corr |>
  as_tibble(
    rownames = "gene_id"
  ) |>
  pivot_longer(
    cols = -gene_id,
    names_to = "Group",
    values_to = "rlog"
  ) |>
  separate(
    col = Group,
    into = c("Treatment", "Replicate")
  ) |>
  mutate(across(Treatment, str_to_upper)) |>
  group_by(gene_id, Treatment) |>
  summarise(
    mean_rlog = mean(rlog),
    .groups = "drop"
  ) |>
  pivot_wider(
    names_from = Treatment,
    values_from = mean_rlog
  ) |>
  select(-gene_id) |>
  ggpairs(
    lower = list(
      continuous = wrap("points", size = 0.1, alpha = 0.5)
    ),
    xlab = "Average expression (rlog)"
  ) +
  theme_few()

```

## Differential expression analysis

Because the batch effect can be associated with the replicate to some extent, we can model it in our analysis using:

$\sim Replicate + Treatment$

Then we ignore the replicate factor to work with batch-modelled treatment effects.

```{r}
#| label: full models
#| message: false

# Model fitting and QC with batch modelling
# 
dds_darpin_true <- DESeqDataSetFromMatrix(
  countData = matrix_darpin_counts,
  colData = table_darpin_colData,
  design = ~ Replicate + Treatment
) |>
  DESeq()

plotDispEsts(dds_darpin_true)

```

Dispersion estimation looks good. Now we can call results.

```{r}
#| label: call results
#| message: false

# Result calling
# 
treatments_darpin <- c("a2", "b1", "d4")
names(treatments_darpin) <- str_to_upper(treatments_darpin)

table_results_darpin <- map_dfr(
  .x = treatments_darpin,
  .f = \(x) get_results(
    dds_object = dds_darpin_true,
    group = x,
    control = "e35",
    table_annotation = table_annotation,
    alpha = 0.01,
    fc_cutoff = 1.5
  ),
  .id = "Treatment"
)

```

For the DARPin experiment, criteria for significant differential expression are:

-   FDR 0.01
-   Fold change more/less than 1.5

This yields so many genes:

```{r}
#| label: DEG counts

table_results_hitcounts <- table_results_darpin |>
  mutate(across(Regulation, \(x) replace_na(x, "Stable"))) |>
  count(Treatment, Regulation) |>
  filter(
    !is.na(Regulation)
  ) |>
  pivot_wider(
    names_from = Regulation,
    values_from = n
  ) |>
  mutate(
    Up_percent = Up / (Up + Down + Stable) * 100,
    Down_percent = Down / (Up + Down + Stable) * 100,
    Up = glue("{Up} ({round(Up_percent, 0)}%)"),
    Down = glue("{Down} ({round(Down_percent, 0)}%)")
  ) |>
  select(Treatment, Up, Down) |>
  mutate(across(everything(), as.character)) |>
  as.data.frame()

table_results_hitcounts

```

## MA plots

MA plots show average expression on the x-axis vs differential expression on the y-axis. They are useful to distinguish between differentially expressed genes with low and high average expression. Here, average expression adds some biological meaning to the results.

```{r}
#| label: MA plots DARPin
#| fig-height: 3
#| fig-width: 7

# MA plots for DARPin experiments
#
table_results_darpin |>
  ggplot(aes(
    x = log10(baseMean),
    y = log2FoldChange,
    color = Significant
  )) +
  geom_point(
    size = 0.3
  ) +
  scale_x_continuous(
    name = "Average expression (log10 mean normalized counts)",
    limits = c(0, 8)
  ) +
  scale_y_continuous(
    name = "Fold change (log2)",
    limits = c(-11, 11)
  ) +
  scale_color_manual(
    values = c(
      `TRUE` = hsv(h = 118/255, s = 145/255, v = 119/255),
      `FALSE` = hsv(h = 160/255, s = 0, v = 181/255)
    )
  ) +
  facet_wrap(~Treatment, nrow = 1) +
  theme_few() +
  theme(
    legend.position = "none"
  )

```

## Pairwise similarity of treatments for DARPin experiment

We visualise the correlation between 3 treatments within the DARPin experiment. This is to compare overall similarity between effects of different DARPin treatments.

```{r}
#| label: DARPin pairwise similarity

# Pairwise similarity for DARPin samples
#
matrix_darpin_lfc_cor <- table_results_darpin |>
  select(Treatment, ensembl_gene_id_version, log2FoldChange) |>
  pivot_wider(
    names_from = Treatment,
    values_from = log2FoldChange
  ) |>
  column_to_rownames("ensembl_gene_id_version") |>
  as.matrix() |> 
  cor()

plot_heatmap_correlation <- Heatmap(
  matrix = matrix_darpin_lfc_cor,
  col = c("gray50", "gray10"),
  name = "Correlation",
  cell_fun = \(j, i, x, y, w, h, c) grid.text(
    label = sprintf("%.2f", matrix_darpin_lfc_cor[i, j]),
    x = x,
    y = y,
    gp = gpar(
      fontsize = 9,
      col = "white"
    )
  ),
  cluster_rows = FALSE,
  cluster_columns = FALSE,
  width = unit(3, "cm"),
  height = unit(3, "cm")
) |> 
  draw()

pdf(
  file = here("plots/plot_heatmap_rnaseq_correlation.pdf"),
  height = 2,
  width = 2.5
)

plot_heatmap_correlation

dev.off()

```

## UpSet plot

```{r}
#| label: upset plots

# Upset plots
# 
list_deg <- map2(
  .x = rep(c("Up", "Down"), each = 3),
  .y = rep(names(treatments_darpin), 2),
  .f = \(x, y) unique(pull(filter(table_results_darpin, Regulation == x & Treatment == y), ensembl_gene_id_version))
)

names(list_deg) <- paste(rep(names(treatments_darpin), 2), rep(c("Up", "Down"), each = 3))

pdf(
  file = here("plots/plot_upset_rnaseq_all.pdf"),
  width = 7,
  height = 5
)

upset(
  data = fromList(list_deg),
  sets = rev(names(list_deg)),
  keep.order = TRUE,
  order.by = "freq",
  set_size.show = TRUE,
  set_size.scale_max = 3500
)

dev.off()

```

## Selection of commonly deregulated genes

Genes that are deregulated in the same direction in all 3 DARPin treatments are marked as Common.

```{r}
#| label: commonly deregulated genes

# Selection of commonly deregulated genes
# 
commonly_deregulated_genes <- table_results_darpin |> 
  filter(
    !is.na(Regulation)
  ) |> 
  filter(
    n() == 3,
    .by = c(ensembl_gene_id_version, Regulation)
  ) |> 
  pull(ensembl_gene_id_version) |> 
  unique() |> 
  sort()

table_results_darpin$Common <- table_results_darpin$ensembl_gene_id_version %in% commonly_deregulated_genes

```

## Save results

```{r}
#| label: save

write_tsv(
  x = table_results_darpin,
  file = here("data/table_results_rnaseq.tsv")
)
```

## Session info

```{r}
#| label: sessionInfo

sessionInfo()

```