---
title: "Exercise: Analysis of a real dataset"
author: "Valerio Licursi"
date: "June 26, 2024"
institute: "IBPM-CNR"
format: html
execute:
  echo: true
---

#### Objective:
Perform statistical tests using the `ggpubr` package.

#### Prerequisites:
- Basic knowledge of R and data frames
- `ggpubr` and `tidyverse` packages installed (`install.packages("ggpubr")` and `install.packages("tidyverse")`)

### Exercise Steps:

#### 1. Install and Load Necessary Packages

```{r}
library(tidyverse)
library(ggpubr)
```

#### 1. Import the dataset

```{r}
gene_exp <- readxl::read_excel("../files_day2/gene_expression_irradiated.xlsx", sheet = 1)
head(gene_exp)
```

The first column doesn't have a name therefore we can `rename` it:

```{r}
gene_exp <- readxl::read_excel("../files_day2/gene_expression_irradiated.xlsx", sheet = 1) %>% 
  rename("sample" = `...1`)
head(gene_exp)
```

#### 3. Data Preparation

```{r}
gene_exp <- pivot_longer(gene_exp, cols = starts_with("rep"), names_to = "rep")
head(gene_exp)

```

#### 4. Create a summary of the data

```{r}
gene_exp_summary <- gene_exp %>%
  group_by(Dose, Time) %>%
  dplyr::summarise( 
    n=n(),
    mean=mean(value, na.rm = TRUE),
    sd=sd(value, na.rm = TRUE)
  ) %>%
  mutate( se=sd/sqrt(n))  %>%
  mutate( ypos = mean+sd+30)

gene_exp_summary
```

#### 5. Visualize the Data with a line plot

```{r}

# Reorder factors
gene_exp$Dose <- ordered(gene_exp$Dose, levels=c("Ctrl", "Gy10"))
gene_exp$Time <- ordered(gene_exp$Time, levels=c("2h", "6h", "24h"))
gene_exp_summary$Dose <- ordered(gene_exp_summary$Dose, levels=c("Ctrl", "Gy10"))
gene_exp_summary$Time <- ordered(gene_exp_summary$Time, levels=c("2h", "6h", "24h"))

# Line plot ----
pd <- position_dodge(0) # avoid errorbar overlap
figure <-
  ggplot(data=gene_exp_summary, aes(x=Time, y=mean, colour=Dose, group=Dose)) + 
  geom_line(size=1.1, position=pd) +
  geom_point(size=4, aes(shape=Dose, fill=Dose), position = pd) + # or +geom_point(size=5, shape=21, fill='white') 
  geom_errorbar(aes(ymin=mean-sd, ymax=mean+sd, colour=Dose), width=0.1, position = pd) + 
  ylab("Intensity Mean (Arbitrary units)") +
  scale_shape_manual(values=c(21,23), labels = c("Ctrl", "10 Gy")) +
  scale_color_manual(values=c("grey40",  "#E41A1C"), labels = c("Ctrl", "10 Gy")) +
  scale_fill_manual(values=c("grey40", "#E41A1C"), labels = c("Ctrl", "10 Gy")) +
  annotate("text", x = c(1,2), y = c(gene_exp_summary$ypos[5], gene_exp_summary$ypos[6]), 
           label = c("****", "***"), 
           size = 6) +
  theme_pubr(base_size = 16, margin = FALSE, legend = "top")
# scale_color_brewer(palette='') per cambiare colori su pattern preimpostati
figure

```

#### 6. Perform ANOVA and post-hoc test

```{r}
TukeyHSD(aov(value ~ sample, data=gene_exp)) %>% broom::tidy() %>% filter(adj.p.value <= 0.05)

```

