9_Localisation_of_PURA_targets.Rmd

---
title: "9 Comparison of PURA targets RNAs enriched in dendrites or cytoplasmic granules"
author: "Melina"
date: "`r format(Sys.time(), '%d %B, %Y')`"
output:
  BiocStyle::html_document:
    toc_float: TRUE
    toc: TRUE0
---

```{r setup, include=FALSE}
require("knitr")
knitr::opts_chunk$set(warning=FALSE, message=FALSE, cache=TRUE) #, fig.pos = "!H", out.extra = ""
```


```{r libraries, include=FALSE}
library(GenomicRanges)
library(knitr)
library(dplyr)
library(ggpubr)
library(biomaRt)
library(purrr)

source("/Users/melinaklostermann/Documents/projects/PURA/02_R_new_pip/XX-helpful-chunks/theme_paper.R")
theme_set(theme_paper())
outpath <- "/Users/melinaklostermann/Documents/projects/PURA/Molitor-et-al-2022/"

```



# What was done?

- Overlap of PURA bound to dendritic RNAs from 6 indipendent RNAseqs of mice dendrites
- 

# Input

```{r}
# binding sites
BS <- readRDS("/Users/melinaklostermann/Documents/projects/PURA/Molitor-et-al-2022/binding_sites_characterized.rds")

# deseq analysis from PURA knockdown
res_RNAseq <- readRDS("/Users/melinaklostermann/Documents/projects/PURA/Molitor-et-al-2022/deseq.rds")
HeLa_expression <- filter(res_RNAseq, baseMean > 0)
rnaseq_counts <- readRDS("/Users/melinaklostermann/Documents/projects/PURA/Molitor-et-al-2022/rnaseq_counts.rds")

# dendritic RNAs
dendritic_RNAs <- readxl::read_xlsx("/Users/melinaklostermann/Documents/projects/PURA/02_R_new_pip/07-Localisation/02-DendriticRNAs/12915_2019_630_MOESM6_ESM.xlsx", sheet = 1)

# p-body transcriptome
pb_transcriptome <- readxl::read_xlsx("/Users/melinaklostermann/Documents/projects/PublicData/p-body/1-s2.0-S1097276517306512-mmc3.xlsx")
# stress granule transcriptome
sg_transcriptome <- readxl::read_xlsx("/Users/melinaklostermann/Documents/projects/PublicData/Stress_granule/1-s2.0-S1097276517307906-mmc2.xlsx")

```

```{r}
# helper functions
give.n <- function(x){
  return(c(y = max(x)*1.1, label = length(x))) 
  # experiment with the multiplier to find the perfect position
}

gg_n <-   stat_summary(fun.data = give.n, geom = "text", fun.y = max, 
                   position = position_dodge(width = 0.75))



```


# Dendritic RNAs
## Studies identifying dendritic RNAs

- There are several studies identifying dendritic RNAs, by performing RNAseq on dendrites (in mice!)
- The study of Middleton et. al 2019 (doi:10.1186/s12915-019-0630-z) offers a list of dendritic RNAs found in 6 different studies.
- We can use genes discribed in 3 or 4 studies as a strong set of dendritic RNAs.
- Here I overlap iCLIP binding sites and RNAseq changing RNAs with dendritic RNAs to suggest some whos dendritic localisation might be influenced by PURA (as a transport factor). Some can be tested for their localisation in neuronal precursors in fish experiments with and without a PURA knockdown.


# How many genes can be compared between mouse dendrites and HeLa cells?

Not all human and mouse genes have one specific otrolog in the other organism. here we only compare those genes that do. Therefore we the mmusculus_homolog_orthology_type = "ortholog_one2one" retrieved from biomart.


```{r}
############
# translate mouse genes to human genes and get ortology information 
#############
ensembl38p13 <- useMart(biomart='ENSEMBL_MART_ENSEMBL',
                        dataset='hsapiens_gene_ensembl')

my_attributes <- c("ensembl_gene_id", "mmusculus_homolog_ensembl_gene", "mmusculus_homolog_associated_gene_name", "mmusculus_homolog_orthology_type", "external_gene_name")

info_biomart_new <- getBM(attributes=my_attributes,

                        mart = ensembl38p13, useCache = F)

# add biomart infos and save files
dendritic_RNAs <- left_join(dendritic_RNAs, info_biomart_new, by = c(GeneName = "mmusculus_homolog_associated_gene_name")) 
save(dendritic_RNAs, file = paste0(outpath, "dendritic_RNAs_with_human_ortholog_one2one.RData"))

HeLa_expression <- left_join(HeLa_expression, info_biomart_new, by = c(gene_id = "ensembl_gene_id")) 
save(HeLa_expression, file = paste0(outpath, "HeLA_RNAs_with_mouse_ortholog_one2one.RData"))

```

Moreover, we only consider genes that were expressed in HeLa (detected in our RNAseq).

```{r}
# add expression
dendritic_RNAs <- mutate(dendritic_RNAs, HeLa_expressed = case_when(
  (ensembl_gene_id %in% HeLa_expression$ensembl_gene_id) & (mmusculus_homolog_orthology_type == "ortholog_one2one") ~ "HeLa_expressed",
                                                                     T ~ "not_HeLa_expressed"))

###################
# plot dendritic perspective
##################

dend_df <- data.frame( all_dend = rep(T,nrow(dendritic_RNAs))) %>% # all dendritic genes
  mutate(human_ortholog = # with a 1-1 orthlog
           case_when(dendritic_RNAs$mmusculus_homolog_orthology_type == "ortholog_one2one" ~ T, T ~ F ),
         HeLa_expression = # expressed in HeLa
           case_when(dendritic_RNAs$HeLa_expressed == "HeLa_expressed" ~ T, T~ F )) 

# make plotable
dend_df <- colSums(dend_df) %>% reshape2::melt(.)
dend_df$variable <- factor(rownames(dend_df), levels =  rownames(dend_df))        

dend_df$value <- c(dend_df$value[1] - dend_df$value[2], dend_df$value[2] - dend_df$value[3], dend_df$value[3] )

ggplot(dend_df)+
  geom_col(aes(x = 1, y = value, fill = variable, position = "stack"))+
   scale_fill_grey(start=0.8, end=0.3)+
  theme(legend.position = "none")

ggsave(paste0(outpath, "dendritic_to_HeLa_RNAs.pdf"), width = 3, height = 6, units = "cm")


###################
# HeLa perspective
##################


hela_df <- data.frame( all_HeLa= rep(T,nrow(HeLa_expression))) %>%
  mutate(Mouse_ortholog = case_when(HeLa_expression$mmusculus_homolog_orthology_type == "ortholog_one2one" ~ T, T ~ F ),
          Dendritic_expression = case_when(HeLa_expression$ensembl_gene_id %in% dendritic_RNAs$ensembl_gene_id ~ T, T~ F )) 

# make plotable
hela_df <- colSums(hela_df) %>% reshape2::melt(.)
hela_df$variable <- factor(rownames(hela_df), levels =  rownames(hela_df))   

hela_df$value <- c(hela_df$value[1] - hela_df$value[2], hela_df$value[2] - hela_df$value[3], hela_df$value[3] )

ggplot(hela_df)+
  geom_col(aes(x = 1, y = value, fill = variable), position = "stack")+
   scale_fill_grey(start=0.8, end=0.3)+
  theme(legend.position = "none")

ggsave(paste0(outpath, "HeLa_to_dendritic_RNAs.pdf"), width = 3, height = 6, units = "cm")

```

## PURA binding in dendritic cells

```{r}

# add PURA binding info to dendritic genes
dendritic_RNAs <- mutate(dendritic_RNAs, bound = case_when(
  ensembl_gene_id %in% BS$gene_id ~ "PURA-bound",
  T ~ "No binding"
) )


##################
# plot num studies perspective
###################

# information of dendritic detection in number of studies 
# on basis on genes recovered in HeLa cells with 1-1- orthologs

dendritic_recovered_in_hela_num_studies <- dendritic_RNAs %>% filter((HeLa_expressed == "HeLa_expressed") & (mmusculus_homolog_orthology_type == "ortholog_one2one"))

ggplot(dendritic_recovered_in_hela_num_studies, aes(x = NumStudiesFound, fill = bound))+
  geom_bar(position = "dodge")+
  scale_y_log10()+
  theme(legend.position = "none")+
  theme_paper()

ggsave(paste0(outpath, "dendritic_binding_HeLa.pdf"), width = 6, height = 6, units = "cm")
  
########################
# plot pie visualisation
########################

# get reproducible found dendritic RNAs
repro_cut  = 3
repro_dendritic_RNAs <- dendritic_RNAs %>% filter((NumStudiesFound >= repro_cut) & mmusculus_homolog_orthology_type == "ortholog_one2one")

# calculate overlaps
venn_df <- data.frame( bound_dend_HeLa = repro_dendritic_RNAs$bound == "PURA-bound",
                       repro_dendritic_RNAs = T)

venn_fit <- eulerr::euler(venn_df)


# make pie chart
gg_df <- venn_fit$original.values %>% reshape2::melt()
gg_df$group = rownames(gg_df)

ggplot(gg_df, aes( x ="", y=value, fill = group))+
  geom_bar(stat="identity", width=1) +
  coord_polar("y", start=0)+
  theme_paper()+
  theme(legend.position =  "None")

ggsave(paste0(outpath, "dendritic_pie.pdf"), width = 5, height = 5, units = "cm")



```

# Cytoplasmic granules

- pb: p-bodies
- sg: stress granules
- the studies from Hubstenberger et al and Khong et al performed RNAseq of pbs or sg vs whole cell extracts and then calculate foldChanges of granule vs whole cell
- we treat RNAs with a foldchange > 0 and a p-value < 0.01 from these studies as enricht in the granule 
- granule enriched RNA is then overlapped with PURA-binding and compared to PURA kd foldchanges

## Overlap of pb and sg enriched RNAs with PURA bound RNAs

```{r}
# rename tables
pb_transcriptome <- mutate(pb_transcriptome, l2fc = as.numeric(`RNA enrichment in P-body (log2) (Fold change=sorted P-bodies/pre-sorted fraction)`), pfdr = as.numeric(`adjusted p-value (FDR) of RNA enrichment`)  )
sg_transcriptome <- sg_transcriptome %>% mutate(l2fc = log2(`Fold change`))


# get pb_transcriptome and sg_transcriptome tpms from our HeLa RNAseq
tpms <- rnaseq_counts[, c("gene_id", "mean_tpm_rnaseq")]
tpms <- filter(tpms, mean_tpm_rnaseq > 0)

pb_transcriptome_expr <- pb_transcriptome[pb_transcriptome$`Ensembl Gene ID` %in% tpms$gene_id,]
sg_transcriptome_expr <- sg_transcriptome[sg_transcriptome$test_id %in% tpms$gene_id,]

pb_transcriptome_expr <- pb_transcriptome_expr %>% left_join(tpms, by= c(`Ensembl Gene ID` = "gene_id"))
sg_transcriptome_expr <- sg_transcriptome_expr %>% left_join(tpms, by= c(test_id = "gene_id"))

# enriched granule RNAs
pb_up_transcriptome <- pb_transcriptome_expr[pb_transcriptome_expr$l2fc > 0 & pb_transcriptome_expr$pfdr < 0.01,]
sg_up_transcriptome <- sg_transcriptome_expr[sg_transcriptome_expr$l2fc > 0 & sg_transcriptome_expr$p_value < 0.01, ]

# bound RNAs enriched in granules
bound_pb_transcriptome <- pb_up_transcriptome[pb_up_transcriptome$`Ensembl Gene ID` %in% substr(BS$gene_id,1,15),]
bound_sg_transcriptome <- sg_up_transcriptome[sg_up_transcriptome$test_id %in% substr(BS$gene_id,1,15),]

# RNAs enriched in both granule types
bound_sg_pb <- bound_pb_transcriptome[bound_pb_transcriptome$`Ensembl Gene ID` %in% bound_sg_transcriptome$test_id, ]
sg_pb<-  pb_up_transcriptome[pb_up_transcriptome$`Ensembl Gene ID` %in% sg_up_transcriptome$test_id, ]
  
# venn of overlaps
all_genes <- c(pb_up_transcriptome$`Ensembl Gene ID`, sg_up_transcriptome$test_id, BS$gene_id) %>% substr(., 1, 15) %>% unique(.)
overlaps2 <- data.frame( bound = all_genes %in% substr(BS$gene_id,1,15),
                         pb = all_genes %in% pb_up_transcriptome$`Ensembl Gene ID`,
                         sg = all_genes %in% sg_up_transcriptome$test_id)
              
 
venn <- eulerr::euler(overlaps2)
plot(venn, quantities = TRUE, shape = "ellipse", fontface = 1)

```

## Compare PURA-binding and p-body/stress granule enrichment 

###  In tpm bins

```{r}
#######################
# pb overlap per tpm
#######################
# tpm distribution comp 
pb_up_transcriptome <- mutate(pb_up_transcriptome, bound = case_when(
  `Ensembl Gene ID` %in% substr(BS$gene_id,1,15) ~ T,
  T~F
))

s <- list(0,1,10,100,1000)

pb_up_transcriptome_tpm_cutoffs <- map(s, ~pb_up_transcriptome %>% filter(mean_tpm_rnaseq > .x))

pb_up_transcriptome_tpm_cutoffs_bound <- map2_dfr(pb_up_transcriptome_tpm_cutoffs, s, ~c(bound = nrow(.x[.x$`Ensembl Gene ID` %in% substr(BS$gene_id,1,15),]),
                                                  all_over_tpm_cut =   nrow(.x),
                                                  tpm_cut = .y))


pb_tpm_cutoffs_bound <- mutate(pb_up_transcriptome_tpm_cutoffs_bound, ratio_bound = bound / all_over_tpm_cut)


ggplot(pb_tpm_cutoffs_bound, aes(x = as.factor(tpm_cut), y = ratio_bound, label = all_over_tpm_cut))+
  geom_col(position = "identity")+
  geom_text(aes(label=paste(round(ratio_bound,2)*100, "% \n of", all_over_tpm_cut)), 
             geom="text", color= "black", size=4)+
  #ggrepel::geom_text_repel()+
  ylim(0,1)+
  theme_paper()

ggsave(paste0(outpath, "pb_and_binding_vs_tpm.pdf"), width = 5 , height  = 5, units = "cm")


#######################
# sg overlap per tpm
#######################
# tpm distribution comp 
sg_transcriptome_expr <- mutate(sg_transcriptome_expr, bound = case_when(
  test_id %in% substr(BS$gene_id,1,15) ~ T,
  T~F
))


sg_transcriptome_up <- sg_transcriptome_expr %>% filter(`Fold change`>0 & q_value < 0.01)

s <- list(0,1,10,100,1000)

sg_transcriptome_up_tpm_cutoffs <- map(s, ~sg_transcriptome_up %>% filter(mean_tpm_rnaseq > .x))

sg_transcriptome_up_tpm_cutoffs_bound <- map2_dfr(sg_transcriptome_up_tpm_cutoffs, s, ~c(bound = nrow(.x[.x$test_id %in% substr(BS$gene_id,1,15),]),
                                                  all_over_tpm_cut =   nrow(.x),
                                                  tpm_cut = .y))



sg_tpm_cutoffs_bound <- mutate(sg_transcriptome_up_tpm_cutoffs_bound, ratio_bound = bound / all_over_tpm_cut)


ggplot(sg_tpm_cutoffs_bound, aes(x = as.factor(tpm_cut), y = ratio_bound, label = all_over_tpm_cut))+
  geom_col(position = "identity")+
  geom_text(aes(label=paste(round(ratio_bound,2)*100, "% \n of", all_over_tpm_cut)), 
             geom="text", color= "black", size=4)+
  #ggrepel::geom_text_repel()+
  ylim(0,1)+
  theme_paper()

ggsave(paste0(outpath, "sg_and_binding_vs_tpm.pdf"), width = 5 , height  = 5, units = "cm")

```

### adjust control group for expression

```{r}
########################
# pb expression controlled negative set
#######################
# introduce tpm bins
pb_transcriptome_expr <- mutate(pb_transcriptome_expr, 
                                bound = `Ensembl Gene ID` %in% BS$gene_id,
                                tpm_bin = cut(log10(mean_tpm_rnaseq), 
                   breaks=c(-Inf, seq(-2, 3.5, 0.25), Inf), 
                   labels=c(1:24)))

ggplot(pb_transcriptome_expr, aes(x = tpm_bin, fill = bound))+
  geom_bar(position = "dodge")+
  ggtitle("tpms bins for pbody RNAs before control matching")

# get distribution of bound genes over tpm bins
bound_bin_distr <- table(pb_transcriptome_expr[pb_transcriptome_expr$bound == T, ]$tpm_bin) %>% as.list()

# get the same amount of random non bound genes per tpm bin
set.seed(42)
unbound_binned_list <- pb_transcriptome_expr[pb_transcriptome_expr$bound == F, ] %>% group_split(tpm_bin, .keep = T)
unbound_matched <- map2(unbound_binned_list, bound_bin_distr, ~.x[sample(1:nrow(.x), size = min(c(nrow(.x), .y))),]) %>% map_dfr(~.x)

# fill up unbound to same n with genes from lower tpm bins
n_bound <- sum(unlist(bound_bin_distr))
n_unbound <- nrow(unbound_matched)
n_missing <- n_bound - n_unbound
n_missing

# get genes not in unbound matched set
unbound_left <- pb_transcriptome_expr[pb_transcriptome_expr$bound == F, ] %>% .[!(.$`Ensembl Gene ID` %in% unbound_matched$`Ensembl Gene ID` ),] 
unbound_fill <- arrange(unbound_left, desc(mean_tpm_rnaseq)) %>% .[1:n_missing,]

pb_transcr_bound_unbound_matched <- rbind(pb_transcriptome_expr[pb_transcriptome_expr$bound == T, ], unbound_matched, unbound_fill )

ggplot(pb_transcr_bound_unbound_matched, aes(x = tpm_bin, fill = bound))+
  geom_bar(position = "dodge")+
  #ggtitle("tpms bins for pbody RNAs after control matching")+
  scale_fill_manual(values =  c("grey", "darkred"))+
  #ggtitle("pbody enrichment in relation to BS site", subtitle = "NA = matched control")+
  theme_paper()+
   theme(legend.position = "None")

########################
# sg expression controlled negative set
#######################
# introduce tpm bins
sg_transcriptome_expr <- mutate(sg_transcriptome_expr,
                                bound = test_id %in% BS$gene_id,
                                  tpm_bin = cut(log10(mean_tpm_rnaseq), 
                   breaks=c(-Inf, seq(-2, 3.5, 0.25), Inf), 
                   labels=c(1:24)))

ggplot(sg_transcriptome_expr, aes(x = tpm_bin, fill = bound))+
  geom_bar(position = "dodge")+
  ggtitle("tpm bins for sg RNAs before control matching")

# get distribution of bound genes over tpm bins
bound_bin_distr <- table(sg_transcriptome_expr[sg_transcriptome_expr$bound == T, ]$tpm_bin) %>% as.list()

# get the same amount of random non bound genes per tpm bin
set.seed(42)
unbound_binned_list <- sg_transcriptome_expr[sg_transcriptome_expr$bound == F, ] %>% group_split(tpm_bin, .keep = T)
unbound_matched <- map2(unbound_binned_list, bound_bin_distr[-1], ~.x[sample(1:nrow(.x), size = min(c(nrow(.x), .y))),]) %>% map_dfr(~.x)

# fill up unbound to same n with genes from lower tpm bins
n_bound <- sum(unlist(bound_bin_distr))
n_unbound <- nrow(unbound_matched)
n_missing <- n_bound - n_unbound
n_missing

# get genes not in unbound matched set
unbound_left <- sg_transcriptome_expr[sg_transcriptome_expr$bound == F, ] %>% .[!(.$test_id %in% unbound_matched$test_id ),] 
unbound_fill <- arrange(unbound_left, desc(mean_tpm_rnaseq)) %>% .[1:n_missing,]

sg_transcr_bound_unbound_matched <- rbind(sg_transcriptome_expr[sg_transcriptome_expr$bound == T, ], unbound_matched, unbound_fill )

ggplot(sg_transcr_bound_unbound_matched, aes(x = tpm_bin, fill = bound))+
  geom_bar(position = "dodge")+
  #ggtitle("tpm bins for sg RNAs before control matching")+
  scale_fill_manual(values =  c("grey", "darkred"))+
  #ggtitle("pbody enrichment in relation to BS site", subtitle = "NA = matched control")+
  theme_paper()+
   theme(legend.position = "None")

```

### Compare enrichment to binding

```{r}
# Comparison pb
ggplot(pb_transcr_bound_unbound_matched, aes(y=l2fc, x = bound, fill = bound))+
  geom_violin()+
  geom_boxplot(width = 0.25)+
  gg_n +
  stat_summary(fun.y=mean, geom="text", show_guide = FALSE, 
               vjust=-0.7, aes( label=round(..y.., digits=1)))+
    stat_compare_means(aes(label = ..p.signif..), label.x = 1.5)+
  scale_fill_manual(values =  c("grey", "darkred"))+
  geom_hline(yintercept = 0)+
  #ggtitle("pbody enrichment in relation to BS site", subtitle = "NA = matched control")+
  theme_paper()+
  theme(legend.position = "None")



# Comparison SG
ggplot(sg_transcr_bound_unbound_matched, aes(y=l2fc, x = bound, fill = bound))+
  geom_violin()+
  geom_boxplot(width = 0.25)+
  gg_n +
  stat_summary(fun.y=mean, geom="text", show_guide = FALSE, 
               vjust=-0.7, aes( label=round(..y.., digits=1)))+
    stat_compare_means(aes(label = ..p.signif..), label.x = 1.5)+
  scale_fill_manual(values =  c("grey", "darkred"))+
  geom_hline(yintercept = 0)+
  #ggtitle("pbody enrichment in relation to BS site", subtitle = "NA = matched control")+
  theme_paper()+
  theme(legend.position = "None")

```


## Comparison of granule localisation to PURA kd RNA changes


```{r}
# get most importanat info from both granules
sg_info <- sg_transcriptome_expr[,c("test_id", "l2fc", "q_value")] 
colnames(sg_info) <- c("gene_id", "l2fc.sg", "fdr.sg")
sg_info_2  <- sg_info %>% filter(l2fc.sg >0 , fdr.sg < 0.01 )

pb_info <- pb_transcriptome_expr[, c("Ensembl Gene ID", "l2fc", "pfdr")]
colnames(pb_info) <- c("gene_id", "l2fc.pb", "fdr.pb")
pb_info_2  <- pb_info %>% filter(l2fc.pb >0 , fdr.pb < 0.01 )

# add PURA kd full
granules <- pb_info %>% 
  left_join(.,res_RNAseq, by = c(gene_id = "ensembl_gene_id")) %>%
  mutate(granule = case_when((l2fc.pb > 0 & fdr.pb < 0.01)  ~ "pb",
                             T ~ "no_granule"))


#plot
ggplot(granules, aes(x = granule, y = log2FoldChange))+
  geom_violin()+
  geom_boxplot(width=0.2)+
   gg_n +
  geom_hline(yintercept = 0)+
  theme_paper()+
    stat_compare_means()
  # ggtitle("PURA kd changes in differnet granules", subtitle = "neg set all RNAs detected sg or pb experiments")

ggsave(paste0(outpath, "violin_pb_vs_not.pdf"), width = 4 , height  = 6, units = "cm")

# add PURA kd full
granules_2 <- full_join(sg_info, pb_info, by = "gene_id") %>% 
  left_join(.,res_RNAseq, by = c(gene_id = "ensembl_gene_id")) %>%
  mutate(granule = case_when((l2fc.pb > 0 & l2fc.sg > 0 & fdr.sg < 0.01 & fdr.pb <0.01)  ~ "sg&pb",
                             (l2fc.pb > 0 & fdr.pb < 0.01)  ~ "pb",
                             (l2fc.sg > 0 & fdr.sg < 0.01 )  ~ "sg",
                             T ~ "no_granule"))

my_comparisons <- list(c("no_granule", "pb"), c("no_granule", "sg"), c("no_granule", "sg&pb"))

#plot
ggplot(granules_2, aes(x = granule, y = log2FoldChange))+
  geom_violin()+
  geom_boxplot(width=0.2)+
   gg_n +
  geom_hline(yintercept = 0)+
  theme_paper()+
    stat_compare_means(comparisons = my_comparisons)+
  ggtitle("PURA kd changes in differnet granules", subtitle = "neg set all RNAs detected sg or pb experiments")

ggsave(paste0(outpath, "violin_pb_vs_not.pdf"), width = 4 , height  = 5, units = "cm")


```