Dissertation_KG_Chapter_1.4.Rmd

---
title: "Role of Piwi-piRNA pathway in somatic and cancer cells"
author: "__Konstantinos Geles__"
date: "Thu Jun 30 2022, Last Update: `r format(Sys.Date(), '%a %b %d %Y')`"
output:
  html_document:
    toc: yes
    toc_depth: 3
    df_print: paged
  pdf_document:
    toc: yes
    toc_depth: 3
  html_notebook: null
editor_options:
  chunk_output_type: console
subtitle: UMG PhD Programme of Molecular and Translational Oncology - Circle XXXIV
---

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

This project contains the scripting part of the Doctoral Dissertation of **Konstantinos Geles** with doi:   

# CHAPTER 1: Role of the PIWI-piRNA pathway in Colorectal Cancer (CRC)  
  
## 1.4 Evaluation of potentially functional piRNA expression in COLO205 

###  piRNA expression in COLO205 cell line after sodium periodate oxidation followed by beta-elimination
For this analysis we have used SPORTS 1.0 to perform all the pre-processing, alignment and quantification steps as shown in the ... GitHub Repository.

I analyze the spike in data to see if the treatment worked and then I search for
the piRNA molecules that have functional modification.

Import Libraries
```{r}
library(readxl)
library(vroom)
library(dplyr)
library(stringr)
library(tidyr)
library(ggplot2)
library(scales)
library(ggpmisc)
library(gridExtra)

```

Import the table with the raw reads
```{r}
piRNA_CRC_cells <- vroom("Chapter_1_4/piRNA_spike_ins_raw.txt") %>% 
    filter(str_detect(smallRNA, "spike")) %>%
    select(smallRNA, starts_with("TO")) %>% 
    pivot_longer(cols = -smallRNA, names_to = "sample", values_to = "read") %>% 
    mutate(smallRNA =  case_when(
               smallRNA == "spikeB1" ~ "SS_22",
               smallRNA == "spikeB2" ~ "SS_28",
               smallRNA == "spikeM3" ~ "mSS_28",
               smallRNA == "spikeM4" ~ "mSS_22"),
           sample = str_replace(sample, "TOTAL", "COLO205")) %>% 
    rename("spike-ins" = smallRNA)
```

PhD theme for plots
```{r}
wes_cols <- c(wesanderson::wes_palettes$BottleRocket2[2:5])

PhD_theme <-
  list(
    scale_fill_manual(values = wes_cols),
    scale_color_manual(values = wes_cols),
    theme_bw() +
      theme(
        panel.border = element_blank(),
        strip.background = element_blank(),
        strip.text = element_text(size = 20, colour = "black"),
        axis.line = element_line(),
        panel.grid.major = element_line(size = 0.2),
        panel.grid.minor = element_line(size = 0.1),
        text = element_text(size = 20),
        axis.title.x = element_text(margin = margin(t = 10, r = 10, b = 10, l = 10),
                                    colour = "black"),
        axis.title.y = element_text(margin = margin(t = 10, r = 10, b = 10, l = 10),
                                    colour = "black"),
        axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 20),
        axis.text.y = element_text(size = 20, colour = "black"),
        plot.title = element_text(hjust = 0.5, colour = "black")
      )
  )
```

plot of spike ins
```{r}
p_treat <- piRNA_CRC_cells %>%
    ggplot() +
    geom_col(mapping = aes(x = sample , y = read, fill = `spike-ins`), position = "fill") +
    ylab("Spike-in Reads") +
    xlab("Samples") +
    ggtitle("Percentage of reads mapped to spike-ins\n for treated and not treated samples of COLO205")+
    scale_y_continuous(labels = scales::percent) +
    PhD_theme

tiff(filename = file.path("FIG_17_spike_in_COLO205.tiff"),
     compression = "none", height = 10, width = 14,  units = 'in', res = 600)
p_treat
dev.off()
```

### check the reads with respect to databases in treated and no treated

import the dataset
```{r}
reads_DBS <- read_xlsx("Chapter_1_3/Table_S3_Sellitto_et_al.xlsx", skip = 4) %>% # remove not annot
    filter(str_detect(Annotation, "TOTAL" )) %>% 
    select(Annotation, miRBase:Unannotated) %>% 
    pivot_longer(cols = -Annotation, names_to = "Database", values_to = "Reads") %>% 
    filter(Database != "rRNAdb")

```

PhD theme for plots
```{r}

PhD_theme <-
  list(
    scale_fill_brewer(palette = "Set1"),
    theme_bw() +
      theme(
        panel.border = element_blank(),
        strip.background = element_blank(),
        strip.text = element_text(size = 20, colour = "black"),
        axis.line = element_line(),
        panel.grid.major = element_line(size = 0.2),
        panel.grid.minor = element_line(size = 0.1),
        text = element_text(size = 20),
        axis.title.x = element_text(margin = margin(t = 10, r = 10, b = 10, l = 10),
                                    colour = "black"),
        axis.title.y = element_text(margin = margin(t = 10, r = 10, b = 10, l = 10),
                                    colour = "black"),
        axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 20),
        axis.text.y = element_text(size = 20, colour = "black"),
        plot.title = element_text(hjust = 0.5, colour = "black")
      )
  )
```

plot of reads after treatment
```{r}
p_DBs_treat <- reads_DBS %>%
    ggplot() +
    geom_col(mapping = aes(x = Annotation , y = Reads, fill = Database), position = "fill") +
    ylab("Mapped Reads") +
    xlab("Sample") +
    ggtitle("Percentage of mapped reads to each sncRNA database \nfor treated and non treated samples of COLO205") +
    scale_y_continuous(labels = scales::percent) +
    PhD_theme

tiff(filename = file.path("FIG_18_treatment_sncRNA_DBs_reads_COLO205.tiff"),
     compression = "none", height = 10, width = 14,  units = 'in', res = 600)
p_DBs_treat
dev.off()
```

### Identifying methylated sncRNAs in COLO205

After running the SPORTS workflow I calculated the normalized with TMM cpm for all the
samples.
We import that table that can be found published in Sellitto et al.
```{r}
sncRNA_treat <- read_xls("Chapter_1_4/Table_S7.xls", skip = 4) %>% # remove not annot
    filter(DB != "noAnnot") %>% 
    mutate(Status = case_when(
        Status == "methylated" ~ "Enriched",
        Status == "not-methylated" ~ "Not Enriched",
        Status == "partially-methylated" ~ "Partly Enriched"
    ))

```

PhD theme for plots
```{r}
wes_cols <- c(wesanderson::wes_palettes$Moonrise2[2],
              "#1C1718",
              wesanderson::wes_palettes$Moonrise2[1]
             )

PhD_theme <-
  list(
    scale_fill_manual(values = wes_cols),
    scale_color_manual(values = wes_cols),
    theme_bw() +
      theme(
        panel.border = element_blank(),
        strip.background = element_blank(),
        strip.text = element_text(size = 20, colour = "black"),
        axis.line = element_line(),
        panel.grid.major = element_line(size = 0.2),
        panel.grid.minor = element_line(size = 0.1),
        text = element_text(size = 20),
        axis.title.x = element_text(margin = margin(t = 10, r = 10, b = 10, l = 10),
                                    colour = "black"),
        axis.title.y = element_text(margin = margin(t = 10, r = 10, b = 10, l = 10),
                                    colour = "black"),
        axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 20),
        axis.text.y = element_text(size = 20, colour = "black"),
        plot.title = element_text(hjust = 0.5, colour = "black")
      )
  )
```

format the table and make main plot
```{r}
number_methylated <-  sncRNA_treat %>%  
    filter( DB != "rRNA") %>% 
    group_by(DB, Status) %>% 
    summarise(Sequences = n()) %>% 
    ungroup()

p_treat_main <- sncRNA_treat %>% 
    filter( DB != "rRNA", Tot_NAIO4_median != 0 ) %>% 
    ggplot(aes(x = log10(Tot_NAIO4_mean + 0.001),
               y = log10(Tot_mean + 0.001), col = Status)) +
    geom_jitter() +
    facet_wrap(facets = "DB", ncol = 2) +
    PhD_theme

tbs <- lapply(split(number_methylated, number_methylated$DB), "[", -1) %>% 
    lapply(FUN = mutate, )

df_plot <- tibble(x = rep(-Inf, length(tbs)), 
             y = rep(Inf, length(tbs)), 
             DB = names(tbs), 
             tbl = tbs)
```

add the tables and make the complete plot
```{r}
p_treat_comp <- p_treat_main + 
    geom_table(data = df_plot, aes(x = x, y = y, label = tbl),
                hjust = 0, vjust = 1, size = 6) +
    ylab("Log10 Mean Counts per Millions, Non-Treated COLO205 Samples") +
    xlab("Log10 Mean Counts per Millions, Treated COLO205 Samples") +
    ggtitle("Effect of sodium-periodate treatment to each of the sncRNA classes") 
    
tiff(filename = file.path("FIG_19_effect_treatment_Colo205_.tiff"),
     compression = "none", height = 12, width = 16,  units = 'in', res = 600)
p_treat_comp
dev.off()
```

### sncRNAs in COLO205 Cytosol and Nucleus

import the datasets for cell fractions and methylated molecules
```{r}
library(edgeR)
dge_contrasts <- readRDS("Chapter_1_4/DBs/2nd_treatment/contrast_list_x1_voomQ_CYT_vs_NUC.rds")
dge_cyt <- vroom("Chapter_1_4/DBs/2nd_treatment/cpm_TMM_CYT_vs_NUC.txt")  

```

format the list to dataframe and filter for methylated molecules

```{r}
dge_contrasts_form_cyt <- dge_contrasts %>% 
    bind_rows(.id = "comparison") %>% 
    filter(comparison %in% c("TreCytvsNuc", "TrevsUntCyt", "UntCytvsNuc")) %>% 
    select(-c(AveExpr, t, P.Value, B)) %>% 
    filter(adj.P.Val< 0.05) %>% 
    select(-adj.P.Val) %>% 
    pivot_wider(names_from = comparison, values_from = c(logFC)) %>% 
    filter(!is.na(TrevsUntCyt)) %>% 
    mutate(log_fc_cyt_nuc = case_when(
        is.na(TreCytvsNuc) & is.na(UntCytvsNuc) ~ 0,
        is.na(TreCytvsNuc) ~ UntCytvsNuc,
        is.na(UntCytvsNuc) ~ TreCytvsNuc,
        TreCytvsNuc > 0 & UntCytvsNuc > 0 ~ UntCytvsNuc,
        TRUE ~ 0
    )) %>% 
    separate(col = smallRNA, into = c("Class", "sncRNA"),sep = "_match_" )  

dge_contrasts_form_nuc <- dge_contrasts %>% 
    bind_rows(.id = "comparison") %>% 
    filter(comparison %in% c("TreCytvsNuc", "TrevsUntNuc", "UntCytvsNuc")) %>% 
    select(-c(AveExpr, t, P.Value, B)) %>% 
    filter(adj.P.Val< 0.05) %>% 
    select(-adj.P.Val) %>% 
    pivot_wider(names_from = comparison, values_from = c(logFC)) %>% 
    filter(!is.na(TrevsUntNuc)) %>% 
    mutate(log_fc_cyt_nuc = case_when(
        is.na(TreCytvsNuc) & is.na(UntCytvsNuc) ~ 0,
        is.na(TreCytvsNuc) ~ UntCytvsNuc,
        is.na(UntCytvsNuc) ~ TreCytvsNuc,
        TreCytvsNuc < 0 & UntCytvsNuc < 0 ~ UntCytvsNuc,
        TRUE ~ 0
    )) %>% 
    mutate(log_fc_cyt_nuc = -(log_fc_cyt_nuc)) %>% 
    separate(col = smallRNA, into = c("Class", "sncRNA"),sep = "_match_" )
```

PhD theme for plots
```{r}
PhD_theme <-
  list(
    scale_fill_brewer(palette = "Set1"),
    scale_color_brewer(palette = "Set1"),
    theme_bw() +
      theme(
        panel.border = element_blank(),
        strip.background = element_blank(),
        strip.text = element_text(size = 20, colour = "black"),
        axis.line = element_line(),
        panel.grid.major = element_line(size = 0.2),
        panel.grid.minor = element_line(size = 0.1),
        text = element_text(size = 20),
        axis.title.x = element_text(margin = margin(t = 10, r = 10, b = 10, l = 10),
                                    colour = "black"),
        axis.title.y = element_text(margin = margin(t = 10, r = 10, b = 10, l = 10),
                                    colour = "black"),
        axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 20),
        axis.text.y = element_text(size = 20, colour = "black"),
        plot.title = element_text(hjust = 0.5, colour = "black")
      )
  )
```

make the plots
```{r}
tbs_cyt <- dge_contrasts_form_cyt %>% 
    filter(!Class %in% c("noAnnot", "rRNA")) %>% 
    group_by(Class) %>% 
    summarise("Methylated Sequences" = n())

plot_cyt <- dge_contrasts_form_cyt %>% 
    filter(!Class %in% c("noAnnot", "rRNA"), log_fc_cyt_nuc >= 0) %>% 
    ggplot() +
    geom_point(mapping = aes(x = log_fc_cyt_nuc , y = TrevsUntCyt, col = Class)) +
    ylab("Non-Treated <----- Log2FC -----> Treated") +
    xlab("Log2FC -----> Cytosol enriched sncRNAs") +
    ggtitle("Effect of sodium-periodate treatment for the Cytosol and Nucleus fractions") +
    annotate(geom = "table", x = 6.5, y = 4, 
             label = list(tbs_cyt), 
             vjust = 0, hjust = 0.5, size = 5) +
    PhD_theme


tbs_nuc <- dge_contrasts_form_nuc %>% 
    filter(!Class %in% c("noAnnot", "rRNA")) %>% 
    group_by(Class) %>% 
    summarise("Methylated Sequences" = n())

plot_nuc <- dge_contrasts_form_nuc %>% 
    filter(!Class %in% c("noAnnot", "rRNA"), log_fc_cyt_nuc >= 0) %>% 
    ggplot() +
    geom_jitter(mapping = aes(x = log_fc_cyt_nuc , y = TrevsUntNuc, col = Class)) +
    ylab("Non-Treated <----- Log2FC -----> Treated") +
    xlab("Log2FC -----> Nucleus enriched sncRNAs") +
    annotate(geom = "table", x = 5.5, y = 6, 
             label = list(tbs_nuc), 
             vjust = 0, hjust = 0.5, size = 5)+
    PhD_theme


tiff(filename = file.path("FIG_20_effect_treatment_and_fractionation_Cyt_Nuc_Colo205_.tiff"),
     compression = "none", height = 12, width = 16,  units = 'in', res = 600)
grid.arrange(plot_cyt, plot_nuc)
dev.off()
```