Dissertation_KG_Chapter_1.1.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.1 Exploration of the PIWIL gene expression data in public datasets 

### *PIWIL* genes in Genotype-Tissue Expression (GTEx) project  
```{r libraries}
library(vroom)
library(dplyr)
library(stringr)
library(tidyr)
library(forcats)
library(ggplot2)
```

Explore GTEx and make FIG 8
Import the data to R
```{r import median gene expression}
gtex_median_v8 <- vroom("Chapter_1_1/GTEx_Analysis_2017-06-05_v8_RNASeQCv1.1.9_gene_median_tpm.gct.gz", skip = 2) %>% 
    filter(str_detect(Description, "PIWIL")) %>% 
    pivot_longer(cols = - c(Name, Description), names_to = "Tissue", values_to = "MEDIAN_TPM") %>% 
    mutate(Tissue = str_replace_all(Tissue, "-", "_") %>% str_replace_all(" ", "")) %>%
    filter(!str_detect(Tissue, "Cells"))

mean_med <- gtex_median_v8 %>% 
    separate(col = Tissue, into = c("Tissue", "subTissue"), sep = "_") %>% 
    group_by(Tissue, Description) %>% 
    summarise(medianTPM = median(MEDIAN_TPM)) %>% 
    ungroup() %>% 
    arrange(desc(medianTPM)) %>% 
    mutate(Tissue = as_factor(Tissue))
```

PhD theme for plots
```{r}
wes_cols <- c(wesanderson::wes_palettes$BottleRocket1)

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")
      )
  )
```

Make the figure 
```{r}
plot_GTEx_PIWIL <- mean_med %>% 
    ggplot(aes(x = Description, y = medianTPM, fill = Description)) +
    geom_col()+
    facet_wrap(facets = "Tissue") +
    xlab("Gene") +
    ylab("Median TPM expression") +
    ggtitle("Bulk Tissue gene expression for PIWIL genes") +
    PhD_theme

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

session info and restart R
```{r}
sessionInfo()
```

>R version 4.1.2 (2021-11-01)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 19043)

>Matrix products: default

>locale:
[1] LC_COLLATE=English_United Kingdom.1252  LC_CTYPE=English_United Kingdom.1252   
[3] LC_MONETARY=English_United Kingdom.1252 LC_NUMERIC=C                           
[5] LC_TIME=English_United Kingdom.1252    

>attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

>other attached packages:
[1] ggplot2_3.3.5 forcats_0.5.1 tidyr_1.2.0   stringr_1.4.0 dplyr_1.0.8   vroom_1.5.7  
 

### *PIWIL* genes in TCGA projects  

Explore recount TCGA 
import libraries
```{r warning=FALSE, include=FALSE}
library(recount3)
library(dplyr)
library(readr)
library(tidyr)
library(furrr)
library(purrr)
library(forcats)
library(edgeR)
library(ggplot2)
library(ggrepel)
library(GGally)
library(ggpubr)
library(ComplexHeatmap)
library(circlize)

library(tidybulk)
library(tidySummarizedExperiment)
```

get TCGA data
```{r}
human_projects <- recount3::available_projects(organism = "human")
project_info <- filter(human_projects, str_detect(file_source,"tcga")) 
projects <- project_info %>% 
    pull(project) %>% 
    set_names()
```

pick per one dataset and make rse object
```{r}
tissues <- map(projects, ~ project_info %>% filter(project == .x) %>% 
                   create_rse(type = "gene", jxn_format = "UNIQUE"))
```

scale and transform counts to TPM
```{r}
for(i in seq_along(tissues)){
assay(tissues[[i]], "counts") <- transform_counts(tissues[[i]])
assays(tissues[[i]])$TPM <- recount::getTPM(tissues[[i]])
}

rm(human_projects)
gc()
tissues %>% write_rds("Chapter_1_1/tissues_TCGA_list_rse.rds")
#tissues <- read_rds("Chapter_1_1/tissues_TCGA_list_rse.rds")
```

pick PIWIL genes names 
```{r}
ensg_names_piwil <- tissues[[1]] %>% 
    rowRanges() %>% 
    filter(str_detect(gene_name,"PIWIL")) %>% 
    as_tibble() %>% 
    select(gene_id, gene_name)
```

pick piwil TPM expression
```{r}
piwil_exp <- lapply(tissues, function(x){assays(x)$TPM %>% 
        as_tibble(rownames = "gene_id") %>% 
        filter(gene_id %in% ensg_names_piwil$gene_id)}) %>% 
    bind_rows(.id = "project") %>% 
    inner_join(ensg_names_piwil, by = c("gene_id")) %>% 
    select(gene_id, gene_name, everything()) 

piwil_exp <- piwil_exp %>% 
    pivot_longer(cols = -c(gene_id, gene_name, project), 
                 names_to = "external_id", values_to = "TPM") %>% 
    drop_na() 

```

coldata
```{r}
metad <- map(tissues, ~colData(.x) %>% 
    as_tibble() %>% 
    select(external_id, project = study, 
           project_name = tcga.gdc_cases.project.name, 
           primary_site = tcga.gdc_cases.project.primary_site,
           condition = tcga.cgc_sample_sample_type,
           platform = tcga.gdc_platform, 
           center = tcga.gdc_center.short_name)) %>% 
    bind_rows()

piwil_exp <- piwil_exp %>% left_join(metad)

piwil_exp %>% vroom::vroom_write("TCGA_PIWIL_expression.txt")
#piwil_exp <- vroom::vroom("Chapter_1_1/TCGA_PIWIL_expression.txt")
```

keep datasets that have: "Primary Tumor" and "Solid Tissue Normal"
```{r}
filtered_projects <- piwil_exp %>% 
    filter(condition %in% c("Primary Tumor", "Solid Tissue Normal")) %>% 
    count(project, condition) %>% 
    count(project) %>% 
    filter(n > 1) %>% 
    pull(project)

filtered_projects # 24 projects

piwil_exp %>% 
    filter(condition %in% c("Primary Tumor", "Solid Tissue Normal"),
           project %in% filtered_projects) %>% 
    count(platform) # 2 platforms

piwil_exp %>% 
    filter(condition %in% c("Primary Tumor", "Solid Tissue Normal"),
           project %in% filtered_projects) %>% 
    count(center) # 3 centers

piwil_exp <- piwil_exp %>% 
    filter(condition %in% c("Primary Tumor", "Solid Tissue Normal"),
           project %in% filtered_projects)

tissues <- tissues[filtered_projects]
```

keep only the columns for DIF EXPR
```{r}
plan(multisession, workers = 2)

tissues <- future_map(tissues, ~.x %>% 
                        tidySummarizedExperiment::select(.feature:study,
                             tcga.gdc_cases.project.name,
                             tcga.gdc_cases.project.primary_site,
                             tcga.cgc_sample_sample_type,
                             tcga.gdc_platform, 
                             tcga.gdc_center.short_name) %>% 
                        tidySummarizedExperiment::select(-raw_counts) %>% 
                        tidySummarizedExperiment::rename(
                             project_name = tcga.gdc_cases.project.name,
                             primary_site = tcga.gdc_cases.project.primary_site,
                             condition = tcga.cgc_sample_sample_type,
                             platform = tcga.gdc_platform, 
                             center = tcga.gdc_center.short_name) %>% 
                        tidySummarizedExperiment::mutate(
                             condition = str_replace_all(condition, " ", "_")) %>% 
                        left_join(ensg_names_piwil, by = c(".feature" = "gene_id")))

# remove non standard samples
tissues <- future_map(tissues, ~.x %>%
                          tidySummarizedExperiment::filter(condition %in% c("Primary_Tumor", 
"Solid_Tissue_Normal")))
## save the list
tissues %>% write_rds("Chapter_1_1/tissues_TCGA_list_rse_for_DE.rds")

## tissues  <- read_rds("Chapter_1_1/tissues_TCGA_list_rse_for_DE.rds")
## Specify design and use tidybulk 
plan(multisession, workers = 2)
explore_batch <- tissues %>%  
    future_map(~.x %>% count(condition, project_name, primary_site, platform, center)) %>% 
    bind_rows(.id = "study") 

# edger
edger_get_DE <- function(x){
   
    
    dgList <- edgeR::SE2DGEList(x)

    group <- factor(dgList$samples$condition)
    keep.exprs <- filterByExpr(dgList, group = group)
    dgList <- dgList[keep.exprs,]
    dgList <- calcNormFactors(dgList)

    design <- model.matrix(~0 + dgList$samples$condition)
    colnames(design) <- colnames(design) %>% 
    str_remove("dgList\\$samples\\$condition") 
    rownames(design) <- rownames(dgList$samples)
    
    con_mat <- makeContrasts( Tumour_Normal = Primary_Tumor - Solid_Tissue_Normal,
                              levels = design)
    
    dgList <- estimateDisp(dgList, design, robust=TRUE)
    dgList <- glmQLFit(dgList, design, robust=TRUE)
    dgList <- glmQLFTest(dgList, contrast = con_mat[,"Tumour_Normal"])
    
    return(dgList)
}

## model matrix with batch
edger_get_DE_batch <- function(x){
   
    
    dgList <- edgeR::SE2DGEList(x)

    group <- factor(dgList$samples$condition)
    keep.exprs <- filterByExpr(dgList, group = group)
    dgList <- dgList[keep.exprs,]
    dgList <- calcNormFactors(dgList)

    design <- model.matrix(~0 + dgList$samples$condition + dgList$samples$platform)
    colnames(design) <- colnames(design) %>% 
    str_remove("dgList\\$samples\\$condition") %>% 
    str_remove("dgList\\$samples\\$platform")
    rownames(design) <- rownames(dgList$samples)
    
    con_mat <- makeContrasts( Tumour_Normal = Primary_Tumor - Solid_Tissue_Normal,
                              levels = design)
    
    dgList <- estimateDisp(dgList, design, robust=TRUE)
    dgList <- glmQLFit(dgList, design, robust=TRUE)
    dgList <- glmQLFTest(dgList, contrast = con_mat[,"Tumour_Normal"])
    
    return(dgList)
}
```

## BLCA, BRCA, CESC, CHOL, ESCA, GBM, HNSC, KICH, KIRC, KIRP, LIHC, LUAD, LUSC, PAAD, PCPG, PRAD, SARC, THCA, THYM
```{r}
no_batch_tissues <- c("BLCA", "BRCA", "CESC", "CHOL", "ESCA", "GBM", "HNSC", 
                      "KICH", "KIRC", "KIRP", "LIHC", "LUAD", "LUSC", "PAAD", 
                      "PCPG", "PRAD", "SARC", "THCA", "SKCM", "THYM")
tictoc::tic()
no_batch_tissues_DE <- tissues[no_batch_tissues] %>% 
    future_map(~edger_get_DE(.x))
tictoc::toc()
```

# COAD, READ, STAD, UCEC
```{r}
batch_tissues <- c("COAD", "READ", "STAD", "UCEC")
tictoc::tic()
batch_tissues_DE <- tissues[batch_tissues] %>% 
    future_map(~mutate(.x, platform = str_replace(platform, " ", "_")) %>% edger_get_DE_batch())
tictoc::toc()
```

# join the datasets
```{r}
tissues <- c(batch_tissues_DE, no_batch_tissues_DE) 
tissues %>% write_rds("Chapter_1_1/tissues_TCGA_dge_list_edger_DE.rds") 
rm(batch_tissues_DE, no_batch_tissues_DE)
gc()
plan(multisession, workers = 1)
```

#### Import the dataset and START from here
```{r}
tissues <- read_rds("Chapter_1_1/tissues_TCGA_dge_list_edger_DE.rds") 
```

check DE per tissue
```{r}
tissues_tidy <- tissues %>% 
    map(~.x %>% topTags(object = ., n = nrow(.), 
           adjust.method = "BH", sort.by = "PValue", p.value = 1) %>%
            .$table %>% 
            as_tibble(rownames = "smallRNA") %>% 
            filter(!is.na(gene_name)) %>% 
    select(smallRNA, gene_name:FDR)) %>% 
    bind_rows(.id = "study")

ensg_names_piwil <- tissues_tidy %>% 
    distinct(smallRNA, gene_name)
```

# get log cpm values for plotting
```{r}
lcpm_piwil <- tissues %>% 
    map(~cpm(.x, log = TRUE, prior.count = 4) %>% 
            as_tibble(rownames = "ENSBLID") %>% 
            filter(ENSBLID %in% ensg_names_piwil$smallRNA) %>% 
            pivot_longer(cols = -ENSBLID, names_to = "external_id", values_to = "lcpm") %>% 
            left_join(as_tibble(.x$samples))) %>% 
    bind_rows() %>% 
    left_join(ensg_names_piwil, by = c("ENSBLID" = "smallRNA"))
```

PhD theme for plots
```{r}
wes_cols <- c(wesanderson::wes_palettes$Darjeeling1, 
              wesanderson::wes_palettes$Rushmore1,
              wesanderson::wes_palettes$IsleofDogs1)

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")
      )
  )
```

pick studies with more than 2 normal samples
```{r}
sel_studies <- lcpm_piwil %>% 
    filter(condition == "Solid_Tissue_Normal") %>% 
    count(study, condition, gene_name) %>% 
    filter(n > 2) %>% 
    count(study) %>% 
    pull(study)
```

make violin plot
```{r}

my_comparisons <- list(c("Primary_Tumor", "Solid_Tissue_Normal"))
# PIWIL1

p <- lcpm_piwil %>%  
    filter(study %in% sel_studies, gene_name == "PIWIL1") %>% 
    mutate(study = as.factor(study), condition = as_factor(condition)) %>%
    ggviolin(x = "condition", y = "lcpm", fill = "study",
             add = "boxplot", 
             add.params = list(fill = "white")) +
    facet_wrap(~study, nrow = 3) + 
    stat_compare_means(comparisons = my_comparisons, label = "p.signif")+
    stat_compare_means(label.y = 10) +
    ggtitle('TCGA PIWIL1 bulk expression difference between\nPrimary Tumor and Solid Tissue samples') +
    xlab("Condition") +
    ylab("log Counts per Million gene expression") +
    PhD_theme

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

```

make heatmap of logFC for PIWIL genes
```{r}
mat_v <- tissues_tidy %>%
    filter(FDR < 0.05) %>% 
    select(study, gene_name, logFC) %>% 
    pivot_wider(names_from = "gene_name", values_from = logFC) %>% 
    tibble::column_to_rownames("study") %>% 
    as.matrix()

sc_mat <- mat_v %>% t() %>% scale() %>% t()
sc_mat %>% dim()
sc_mat %>% head()
hist(sc_mat)
hist(mat_v)

f <- colorRamp2(c(-1, 0 ,1),c("#53868B", "#8B8878", "#FFD700"))

ha_1 <- HeatmapAnnotation(GENES = colnames(mat_v),
            col = list(GENES =
                           wesanderson::wes_palettes$Royal1[1:4] %>% 
                           set_names(colnames(mat_v))))

ht_1 <- Heatmap(
        matrix = mat_v, #data
        top_annotation = ha_1, #annot
        col = f, #colors data  
        na_col = "grey",
        
        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)),
        cluster_columns = FALSE,
        cluster_rows = FALSE,
        show_row_dend = FALSE,
        show_row_names = TRUE,
        show_column_names = TRUE,
        name = "Scaled\nlogFC",
        clustering_distance_columns = "spearman",
        clustering_method_columns =  "ward.D2",
        column_order = c(1,4,3,2), 
        row_title = "Study",
        column_title = "Heatmap of significantly DE\nlog fold changes of PIWIL genes per study"
        
)

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

### *PIWIL* genes in Expression Atlas, E-MTAB-2706 project  

Import Libraries
```{r}
library(vroom)
library(dplyr)
library(stringr)
library(tidyr)
library(forcats)
library(ggplot2)
```

Using the url: we download the metadata regarding each cell line.
```{r}
cell_lines_metdat <- vroom("Chapter_1_1/E-MTAB-2706.sdrf.txt") %>% 
    select(Source_Name = "Source Name", Cell_line = "Characteristics[cell line]",
           tissue = "Characteristics[organism part]",
           disease = "Characteristics[disease]") %>% 
    filter(str_detect(tissue, "colon")) %>% 
    mutate(Source_Name = str_replace(Source_Name, " ","_")) %>% 
    distinct(Source_Name, Cell_line, .keep_all = TRUE)

gene_data <- vroom("Chapter_1_1/E-MTAB-2706.additional.1/140625_Klijn_geneToTranscript.txt") %>% 
    rename(geneID = "entrez")
```

Using the data on the url: we import the expression values for each cell line.
```{r}
RPKM_cell_lines <- vroom("Chapter_1_1/E-MTAB-2706.additional.1/140625_Klijn_RPKM_coding.txt") %>% 
    pivot_longer(cols = -geneID, names_to = "samples", values_to = "RPKM") %>% 
    mutate(Source_Name = str_replace(samples, " ", "_")) %>% 
    inner_join(cell_lines_metdat) %>% 
    inner_join(gene_data)

PIWIL_genes_expression <- RPKM_cell_lines %>% 
    filter(str_detect(gene_name, "PIWIL")) 
```

PhD theme for plots 
```{r}
wes_cols <- c(wesanderson::wes_palettes$Rushmore[1:4],
                             wesanderson::wes_palettes$Darjeeling1)

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")
      )
  )
```

Make the plot
```{r}

p_pw_rpkm_piwils <- PIWIL_genes_expression %>% 
    filter(
           Cell_line %in% c("Caco-2", "COLO 205", "HCT 116", "HT-115",
                            "HT-29", "SW 1417", "SW 403", "RKO")) %>% 
    ggplot() +
    geom_col(mapping = aes(x = Cell_line, y = RPKM, fill = Cell_line))+
    facet_wrap(facets = "gene_name", ncol = 1, scales = "free_y") +
    xlab("Cell Line") +
    ggtitle("Expression of PIWIL genes in CRC cell-lines \n from the Public dataset: E-MTAB-2706") +
    PhD_theme

tiff("FIG_12_E_MTAB_2706_PIWILs_rpkm.tiff", width = 9, height = 8, units = 'in', res = 300)
p_pw_rpkm_piwils
dev.off()
```