Dissertation_KG_Chapter_2.3.1.Rmd

---
title: "Role of Piwi-piRNA pathway in somatic and cancer cells"
author: "__Konstantinos Geles__"
date: "Wed Jul 13 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 2: Data Analysis Workflow for small-RNAseq focused on piRNAs  
  
## 2.3.1 Analysis of samples from patients with CRC and comparison with previous analyses.

Import Libraries
```{r}
library(vroom)           #
library(readr)           #
library(dplyr)           #
library(tidyr)           # data import and manipulation
library(purrr)           #
library(tibble)          #
library(stringr)         #
library(forcats)         #

library(ggplot2)         #
library(scales)          # make exploration data plots
library(ggpmisc)         #

library(edgeR)           # identify DEs and
library(ComplexHeatmap)  # Heatmaps
library(circlize)        #

library(clusterProfiler) #
library(wppi)            # Functional enrichment analysis and 
library(ReactomePA)      # protein protein interaction network
library(enrichplot)      #
```

###  DE sncRNAs per class

import the resulted de sncRNAs of CRC patients
```{r}
CRC_tissues <- vroom("Chapter_2_3/DEA_piRNA_CC_tissues_GRCh38_25_Jan_2022/all_comparisons_long_voom_TMM_salmon_fc_LFCs_25_Jan_2022.txt")

```

PhD theme for plots
```{r}
wes_cols <- c(wesanderson::wes_palettes$Cavalcanti1,
              wesanderson::wes_palettes$Royal2,
              wesanderson::wes_palettes$Darjeeling2)

PhD_theme <- # this PHD theme !!!!
  list(
    scale_fill_manual(values = wes_cols),
    theme_bw() +
      theme(
        panel.border = element_blank(),
        strip.background = element_blank(),
        strip.text = element_text(size = 25, colour = "black",
                                  face = "bold"),
        axis.line = element_line(),
        panel.grid.major = element_line(size = 0.3),
        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", face = "bold"),
        axis.title.y = element_text(margin = margin(t = 10, r = 10, b = 10, l = 10),
                                    colour = "black", face = "bold"),
        axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 20,
                                   colour = "black", face = "bold"),
        axis.text.y = element_text(size = 20, colour = "black", face = "bold"),
        plot.title = element_text(hjust = 0.5, colour = "black", face = "bold"),
        legend.position = "bottom", 
        legend.direction = "horizontal"
      )
  )
```

plot the amount of DE per sncRNA class
```{r}
plot_CRC_DE <- CRC_tissues %>% 
    filter(adj.P.Val < 0.05) %>% 
    ggplot(aes(x = gene_type, fill = gene_type)) +
    geom_bar(stat = "count") +
    stat_count(geom = "text", colour = "black", size = 7.5, 
              aes(label = ..count..), position  = position_fill(vjust = 255 )) +
    facet_wrap(facets = "quantification", ncol = 1) + 
    xlab("sncRNA Class") +
    ylab("Count of different sequences") +
    ggtitle("Amount of DE sncRNAs with each quantification\n approach in CRC patients in-house dataset") +
    PhD_theme

tiff(filename = file.path("FIG_40_CRC_DE_salmon_FC.tiff"),
     compression = "none", height = 12, width = 16,  units = 'in', res = 600)
plot_CRC_DE
dev.off()
```

### DE piRNAs and heatmap of union of both approaches

DE piRNAs, keep the union of both approaches
```{r}
DE_piRNA <- CRC_tissues %>% 
    filter(adj.P.Val < 0.05, gene_type == "piRNA") %>% 
    select(smallRNA:logFC) %>% 
    pivot_wider(names_from = quantification, values_from = logFC) 
```

import the exprs values from featurecounts for the common DE piRNAs
```{r}
exprs_mat_tot <- read_rds("Chapter_2_3/EDA_piRNA_CC_tissues_GRCh38_24_Jan_2022/featureCounts/list_norm_dgls_featureCounts.rds") %>% 
    magrittr::extract2("TMM") %>% #get only the TMM normalized expressions
    edgeR::cpm(log = TRUE, prior.count = 4) %>% # transform to lcpm
    .[DE_piRNA$smallRNA,] # subset to DE piRNAs

```

import the groups table
```{r}
table_groups <- read_rds("Chapter_2_3/EDA_piRNA_CC_tissues_GRCh38_24_Jan_2022/featureCounts/list_norm_dgls_featureCounts.rds") %>% 
    magrittr::extract2("TMM") %>% 
    magrittr::extract2("colours")
```

Scale the matrix for the heatmap
```{r}
# make the exprs matrix for the heatmap -----
sc_mat <- exprs_mat_tot %>% t() %>% scale() %>% t()

sc_mat %>% dim()
sc_mat %>% head()
str_c("min = ", round(min(sc_mat),3),
      "| median = ", round(median(sc_mat),3), "| max = ", round(max(sc_mat),3))

stopifnot(identical(as.character(table_groups$name), colnames(sc_mat)))

# make the logFC matrix for the heatmap -----
scfc_mat <- DE_piRNA %>% column_to_rownames("smallRNA") %>% as.matrix()
scfc_mat %>% dim()
scfc_mat %>% head()
str_c("min = ", round(min(scfc_mat, na.rm = T),3),
      "| median = ", round(median(scfc_mat, na.rm = T),3), "| max = ", round(max(scfc_mat, na.rm = T),3))

scfc_mat <- scfc_mat[rownames(sc_mat),]

rownames(sc_mat) <- rownames(sc_mat) %>% str_remove("_GR.+")

rownames(scfc_mat) <- rownames(scfc_mat) %>% str_remove("_GR.+")

identical(rownames(sc_mat), rownames(scfc_mat))
```

make the heatmap
```{r}
f <- colorRamp2(c(-2, 0, 2), c("#53868B", "#8B8878", "#FFD700"))

f_LFCs <- colorRamp2(c(-2, 0, 2),
                    c("#008B8B", "black", "#B22222"))
# mat
ha_1 <- HeatmapAnnotation(
                          Group = table_groups$group,
                          Status = table_groups$MSI_status,
                          Location = table_groups$Location,
                          BRAF = table_groups$BRAF,
                          KRAS = table_groups$KRAS,
                          Lynch = table_groups$Lynch,
                          MLH1 = table_groups$MLH1,
                          
            col = list(
                Group = table_groups %>%  
                           select(group, group_col) %>% 
                           deframe(),
                Status = table_groups %>%  
                           select(MSI_status, MSI_status_col) %>% 
                           deframe(),
                Location = table_groups %>%  
                           select(Location, Location_col) %>% 
                           deframe(),
                BRAF = table_groups %>%  
                           select(BRAF, BRAF_col) %>% 
                           deframe(),
                KRAS = table_groups %>%  
                           select(KRAS, KRAS_col) %>% 
                           deframe(),
                Lynch = table_groups %>%  
                           select(Lynch, Lynch_col) %>% 
                           deframe(),
                MLH1 = table_groups %>%  
                           select(MLH1, MLH1_col) %>% 
                           deframe() ),
            annotation_name_side = "left",
            annotation_legend_param = 
                list(title_gp = gpar(fontsize = 15, fontface = "bold"),
                     labels_gp = gpar(fontsize = 15, fontface = "bold"))
            )

# lFCS
ha_1_LFCs <- HeatmapAnnotation(
    Method = colnames(scfc_mat),
    col = list(Method = wesanderson::wes_palettes$Moonrise1[c(1,4)] %>%
                   set_names(colnames(scfc_mat))),
    annotation_legend_param = list(title_gp = gpar(fontsize = 15, fontface = "bold"),
                             labels_gp = gpar(fontsize = 15, fontface = "bold")))     

ht_1 <- Heatmap(matrix = sc_mat, #data
        top_annotation = ha_1, #annot
        col = f, #colors data  
        row_title_gp = gpar(fontsize = 20),
        column_title_gp = gpar(fontsize = 20),
        column_names_gp = gpar(fontsize = 20),
        column_names_rot = 45,
        heatmap_legend_param = list(
            title_gp = gpar(fontsize = 15, fontface = "bold"),
            labels_gp = gpar(fontsize = 15)),
        
        show_row_dend = TRUE,
        show_row_names = FALSE,
        show_column_names = FALSE,
        name = "z-score \nabundance",
        column_km = 4,
        row_km = 3,
        clustering_distance_columns = "euclidean",
        clustering_method_columns =  "ward.D2",
        clustering_method_rows = "ward.D2",
        clustering_distance_rows = "euclidean",
        row_dend_reorder = TRUE,
        row_title = "piRNAs",
        column_title = str_glue("Heatmap of {nrow(sc_mat)} common piRNAs\nidentified DE in transcriptomic and genomic approaches"),
)

ht_1_lFCs <- Heatmap(matrix = scfc_mat, #data
        top_annotation = ha_1_LFCs, #annot
        col = f_LFCs, #colors data   
        row_title_gp = gpar(fontsize = 20),
        column_title_gp = gpar(fontsize = 20),
        column_names_gp = gpar(fontsize = 20),
        heatmap_legend_param = list(
            title_gp = gpar(fontsize = 15, fontface = "bold"),
            labels_gp = gpar(fontsize = 15)),
        show_row_dend = FALSE,
        show_column_dend = FALSE,
        show_row_names = FALSE,
        show_column_names = FALSE,
        name = "Log Fold Change",
        clustering_distance_columns = "spearman",
        clustering_method_columns =  "ward.D2",
        clustering_method_rows = "ward.D2",
        clustering_distance_rows = "spearman",
        row_dend_reorder = TRUE
)
draw(ht_1 + ht_1_lFCs,
     column_title = str_glue("Heatmap of {nrow(sc_mat)} DE piRNAs"), merge_legend = TRUE,
     heatmap_legend_side = "bottom", annotation_legend_side = "bottom")

tiff(filename = file.path("FIG_41_piRNA_union_DE_CRC_patients.tiff"),
     compression = "none", height = 10, width = 16,  units = 'in', res = 600)

draw(ht_1 + ht_1_lFCs , merge_legend = TRUE,
     heatmap_legend_side = "bottom", annotation_legend_side = "bottom")

dev.off()

```

### Gene predicted targets for the DE piRNA from CRC patients

import the table with predicted targets
```{r}
total_pred_targ <- vroom("Chapter_2_3/piRNA_predicted_Targets.v02.txt")
```

search the union of DE piRNAs from both FC and salmon
```{r}
union_DE_piRNA <- CRC_tissues %>% 
    filter(adj.P.Val < 0.05, gene_type == "piRNA") %>% 
    distinct(smallRNA) %>% 
    deframe()

total_pred_targ %>% 
    filter(piRNA_id %in% union_DE_piRNA) %>% 
    distinct(piRNA_id)

piRNA_DE_targets <- total_pred_targ %>% 
    filter(piRNA_id %in% union_DE_piRNA) %>% 
    distinct(piRNA_id, Target_gene_name, .keep_all = TRUE) 

piRNA_DE_targets %>% 
    vroom_write("Chapter_2_3/Pred_targets_140_DE_piRNA.txt")

gene_targets <- unique(piRNA_DE_targets$Target_gene_name) 
```

find PPIs for genes of interest
```{r}
# HPO annotations set
HPO_data <- wppi_hpo_data()

# search for gene symbols correlated to cancer of carcinoma
HPO_interest <- HPO_data %>% 
  filter(str_detect(Name, "cancer|carcinoma"))

db <- wppi_data(datasets = c('omnipath', 'kinaseextra')) # here

graph_op <- graph_from_op(db$omnipath)

# subset the graph with regard to the predicted target genes
graph_op_targets <- subgraph_op(graph_op, gene_targets)
igraph::vcount(graph_op_targets)

w_adj_targets <- weighted_adj(graph_op_targets, db$hpo, HPO_interest) 

rw_targets <- random_walk(w_adj_targets)

scores_targets <- prioritization_genes(graph_op_targets, rw_targets, gene_targets)

wppi_scores_targets <- unique(scores_targets$gene_symbol) # 54 genes

total_network_targets <- union(wppi_scores_targets, gene_targets)

```

transform Gene Symbols to ENTREZ IDs
```{r}
Entr_Gene_Ids <- bitr(total_network_targets, 
                     fromType="SYMBOL", 
                     toType="ENTREZID", OrgDb="org.Hs.eg.db")

total_network_targets %>% length()
Entr_Gene_Ids %>% nrow() #75/76
```

over-representation analysis of Ontological databases
```{r}
goenr <- enrichGO(Entr_Gene_Ids$ENTREZID, ont = "BP", OrgDb = 'org.Hs.eg.db', 
                  readable = TRUE) %>%
    as_tibble()

goenr %>% vroom_write("Chapter_2_3/Pred_targets_Enr_GO.txt")

wikienr <- enrichWP(Entr_Gene_Ids$ENTREZID, organism = "Homo sapiens", ) %>%
    as_tibble()

wikienr %>% vroom_write("Chapter_2_3/Pred_targets_Enr_WIKIPATH.txt")

keggenr <- enrichMKEGG(gene = Entr_Gene_Ids$ENTREZID,
                   organism = "hsa",
                   pvalueCutoff = 1,
                   qvalueCutoff = 1) %>%
    as_tibble()

reactome_path <- enrichPathway(gene = Entr_Gene_Ids$ENTREZID, 
                               pvalueCutoff = 0.05, readable = TRUE) %>%
    as_tibble()

reactome_path %>% vroom_write("Chapter_2_3/Pred_targets_Enr_REACTOME.txt")
```

join the enriched pathways
```{r}
wikienr <- wikienr %>% 
    mutate(richFactor = Count / as.numeric(sub("/\\d+", "", BgRatio)))

reactome_path <- reactome_path %>% 
    mutate(richFactor = Count / as.numeric(sub("/\\d+", "", BgRatio)))
```

PhD theme for plots
```{r}
wes_cols <- c(wesanderson::wes_palettes$Zissou1[5:1])
PhD_theme <-
  list(
    scale_fill_brewer(palette = "Set1"),
    theme_bw() +
      theme(
        panel.border = element_blank(),
        strip.background = element_blank(),
        strip.text = element_text(size = 25, colour = "black",
                                  face = "bold"),
        axis.line = element_line(),
        panel.grid.major = element_line(size = 0.3),
        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", face = "bold"),
        axis.title.y = element_text(margin = margin(t = 10, r = 10, b = 10, l = 10),
                                    colour = "black", face = "bold"),
        axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 20,
                                   colour = "black", face = "bold"),
        axis.text.y = element_text(size = 20, colour = "black", face = "bold"),
        plot.title = element_text(hjust = 0.5, colour = "black", face = "bold")
      )
  )
```

plot pathways
```{r}
reactomeenr_plot <- reactome_path %>% 
  head(20) %>% 
  ggplot(
  aes(richFactor, fct_reorder(Description, richFactor))) + 
  geom_segment(aes(xend = 0, yend = Description)) +
  geom_point(aes(color = p.adjust, size = Count)) +
  scale_color_gradientn(colours = wes_cols,
                        guide = guide_colorbar(reverse = TRUE, order = 1)) +
  scale_size_continuous(range=c(2, 10)) +
  #facet_wrap(facets = "Pathway", ncol = 1) +
  xlab("Rich Factor") +
  ylab("Enriched Pathways") + 
  ggtitle("Top 20 Enriched Reactome Pathways in CRC\n patients of the predicted piRNA gene targets")+
  PhD_theme 

tiff(filename = file.path("FIG_42_Enriched_REACTOME_CRC_targets.tiff"),
     compression = "none", height = 12, width = 18,  units = 'in', res = 600)
reactomeenr_plot
dev.off()

wikienr_plot <- wikienr %>% 
  head(20) %>% 
  ggplot(
  aes(richFactor, fct_reorder(Description, richFactor))) + 
  geom_segment(aes(xend = 0, yend = Description)) +
  geom_point(aes(color = p.adjust, size = Count)) +
  scale_color_gradientn(colours = wes_cols,
                        guide = guide_colorbar(reverse = TRUE, order = 1)) +
  scale_size_continuous(range=c(2, 10)) +
  #facet_wrap(facets = "Pathway", ncol = 1) +
  xlab("Rich Factor") +
  ylab("Enriched Pathways") + 
  ggtitle("Top 20 Enriched Wiki Pathways in CRC\n patients of the predicted piRNA gene targets")+
  PhD_theme 


tiff(filename = file.path("FIG_43_Enriched_WIKIPATH_CRC_targets.tiff"),
     compression = "none", height = 12, width = 18,  units = 'in', res = 600)
wikienr_plot
dev.off()
```