3_Characterise_PURA_binding.Rmd

---
title: "3 Characterisation of PURA binding"
author: "Melina"
date: "`r format(Sys.time(), '%d %B, %Y')`"
output:
  BiocStyle::html_document:
    toc_float: TRUE
    toc: TRUE

---

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

# What is done here?

- Types of bound genes
- Distribution of binding sites over gene regions (with and without normalisation for region length)
- Distribution of crosslink event over gene regions, metagene profile 
- Distribution of binding sites per gene and in relation to TPM
- Relative maps of crosslink distribution
- GeneOntology and REACTOME enrichment analysis



```{r libraries, include=FALSE}
library(GenomicRanges)
library(rtracklayer)
library(knitr)
library(GenomicFeatures)
library(dplyr)
library(ggpubr)
library(purrr)
library(wesanderson)
library(cliProfiler)
library(hypeR)


report_color <- (pals::ocean.solar(15))

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

``` 

# Input

```{r input}

# Clip data
binding_sites <-  readRDS("/Users/melinaklostermann/Documents/projects/PURA/Molitor-et-al-2022/binding_sites.rds")  
CL_list <- readRDS("/Users/melinaklostermann/Documents/projects/PURA/Molitor-et-al-2022/bw_merges.rds")



# RNAseq
rnaseq_counts <- readRDS("/Users/melinaklostermann/Documents/projects/PURA/Molitor-et-al-2022/rnaseq_counts.rds")


# annotations
gff3 <- "/Users/melinaklostermann/Documents/projects/anno/GENCODEv31-p12/gencode.v31.annotation.gff3" # for cliProfiler
annotation <- readRDS("/Users/melinaklostermann/Documents/projects/PURA/Molitor-et-al-2022/annotation.rds")
anno_txdb <- loadDb("/Users/melinaklostermann/Documents/projects/PURA/Molitor-et-al-2022/annotation_txdb.db")

```


# Bound Genes and regions

## assign bound gene type 
 
- many binding sites overlap with multiple genes (non codings!)
- in IGV it look like rather the non-coding RNA than the protein-coding RNA is bound in that case
- we therefore use a hieracry on gene types: lncRNA > miRNA > miscRNA > snRNA > snoRNA > protein-coding


 
```{r}
############################
# annote binding sites with the bound ( = overlapping) gene
################################
#  overlaps are resolved by:
# 1. for two genes of different type by the type hieracy
# 2. for to protein coding genes: in this chunk one protein coding gene is assigned at random in the next chunk the region hieracry is used to make this decision and the gene info is changed then

 
annotation_gene <- annotation[annotation$type=="gene"]

# overlaps of bs with any gene
ol <- findOverlaps(annotation_gene, binding_sites, ignore.strand=F)

# index binding sites
binding_sites$idx <-  1:NROW(binding_sites)

# get all possible regions overlapping with crosslinked nucleotide in a temporary file
binding_sites_temp1 <- binding_sites[subjectHits(ol)]
binding_sites_temp1 <- sortSeqlevels(binding_sites_temp1)
elementMetadata(binding_sites_temp1) <-  c(elementMetadata(binding_sites_temp1), elementMetadata(annotation_gene[queryHits(ol), c("gene_id", "gene_type", "gene_name")]))

table(binding_sites_temp1$gene_type)

# chose for each binding site a gene from the highest possible hierarchy
# if two genes from the same type would be present the first one is chosen (randomly)
binding_sites_temp2 <- binding_sites_temp1 %>%  
  as.data.frame(.) %>% # this file contains several entries for 1 binding site overlapping with more than one gene
  group_by(idx) %>% # the idx colum is an index of the binding sites
  arrange(factor(gene_type, levels = c("lncRNA", "miRNA", "miscRNA", "snRNA", "snoRNA", "protein_coding" )), .by_group = T) %>% # arrange by hierarcy
  dplyr::slice(1) # choose randomly  the first gene (the chosen gene is changed in the next chunk for protein coding genes)

binding_sites_temp2 <- makeGRangesFromDataFrame(binding_sites_temp2, keep.extra.columns = T)

table(binding_sites_temp2$gene_type)

# non coding rnas "fake" region non coding
bs_temp_non_cod <- binding_sites_temp2[binding_sites_temp2$gene_type!= "protein_coding",] %>% makeGRangesFromDataFrame(., keep.extra.columns = T)
bs_temp_non_cod$type <- "non_coding"


```

## Binding site region by hiercary approach

- Overlap Bs with mutiple regions and choose one with highest priority
- hieracry: "three_prime_UTR" > "five_prime_UTR" > "CDS" > "intron" 
- only done for protein-coding genes

## Total number of binding sites per region

```{r}

# classify only protein coding genes, regions only from bound genes
anno_regions <- annotation[annotation$gene_type == "protein_coding" & annotation$type %in% c("three_prime_UTR", "five_prime_UTR", "CDS") ]

#get introns
introns <- intronsByTranscript(anno_txdb, use.names=T) %>% unlist(.)
introns$type = "intron"

anno_regions <- c(anno_regions, introns)

#overlaps of binding sites with protein coding genes
binding_sites_temp_pc <- binding_sites_temp2[binding_sites_temp2$gene_type=="protein_coding"]


ol_bs <- findOverlaps(anno_regions, binding_sites_temp_pc-2, ignore.strand = FALSE)


# get all possible regions overlapping with bs
binding_sites_temp_pc2 <- binding_sites_temp_pc[subjectHits(ol_bs)]
binding_sites_temp_pc2 <- sortSeqlevels(binding_sites_temp_pc2)

# get possible region and transcript annotation of bs
anno_bs <- elementMetadata(anno_regions[queryHits(ol_bs)]) %>%
  as.data.frame(.) %>%
  dplyr::select(c("type", "gene_id"))

# and add to temp binding sites
colnames(anno_bs) <- c("type", "gene_id.from_region")
elementMetadata(binding_sites_temp_pc2) <- c(elementMetadata(binding_sites_temp_pc[subjectHits(ol_bs)]), anno_bs)


# choose type highest in hiearcy
binding_sites_temp_pc3 <- as.data.frame(binding_sites_temp_pc2) %>%
    mutate(., gene_id = gene_id.from_region ) %>% # this throws out regions that overlap with non-coding RNAs, with this line no binding sites in introns!!!
  group_by(idx) %>% # specific id per binding site
  arrange(factor(type, levels = c("three_prime_UTR", "five_prime_UTR", "CDS", "intron")), .by_group = T) %>%
  dplyr::slice(1) # arrange by region, select first, transcript info is ignored here

binding_sites_temp_pc3 <- makeGRangesFromDataFrame(binding_sites_temp_pc3, keep.extra.columns = T)
sortSeqlevels(binding_sites_temp_pc3)

# merge BS from non-coding RNAs and binding sites with specific regions
binding_sites_with_regions <- c(binding_sites_temp_pc3, bs_temp_non_cod)


# some protein coding genes not have a region, because no trustworth transcript is annotated in this position of the gene
bs_without_region <- binding_sites[-queryHits(findOverlaps(binding_sites, binding_sites_with_regions))]


```
## Plot total number of gene type

```{r}
###################
# make plots
##################
	
# gene types

bound_genes_gg <- as.data.frame(binding_sites_with_regions) %>%
  filter(!duplicated(gene_id))

ggplot(bound_genes_gg, aes(x = gene_type))+
  geom_bar(fill = "darkgrey")+
  scale_y_log10()+
  coord_flip()+
  theme_paper()

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

# gene regions

binding_sites_gg <- as.data.frame(binding_sites_with_regions)

ggplot(binding_sites_gg, aes(x = type))+
  geom_bar()+
  stat_count(aes(label=paste0(
    sprintf("%1.1f", ..count../sum(..count..)*100),
                              "% \n", ..count..), y=0.5*..count..), 
             geom="text", colour="white", size=4, position=position_dodge(width=1)) +
  theme_paper()+
  ylab("Number of crosslink events")

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

```


```{r}
# final numbers
table(binding_sites_with_regions$region)
table(binding_sites_with_regions[!duplicated(binding_sites_with_regions$gene_id),]$gene_type)

NROW(binding_sites_with_regions[which(!duplicated(binding_sites_with_regions$gene_id)),])

# clean up unnecessary columns
colnames(elementMetadata(binding_sites_with_regions)) <- c(colnames(elementMetadata(binding_sites_with_regions))[1:10], "region", "")
elementMetadata(binding_sites_with_regions) <- elementMetadata(binding_sites_with_regions)[c(1:11)]

```


## Normlisation for region length

```{r}
# classify only protein coding genes
anno_regions_r <- GenomicRanges::reduce(anno_regions)

# annotate the whole annotation to a specific region by hierarcy

# get all 3'utrs
utr3 <- anno_regions[anno_regions$type=="three_prime_UTR"] %>%
  GenomicRanges::reduce()
# remove 3'utr regions from genome coordiantes left to assign
anno_left <- GenomicRanges::setdiff(anno_regions_r, utr3, ignore.strand=FALSE)

# get 5'UTRs from genome coordiantes left to assign
utr5 <- GenomicRanges::intersect(anno_left,
                                 anno_regions[anno_regions$type=="five_prime_UTR"] ) %>%
  GenomicRanges::reduce()
# remove 5'utr regions from genome coordiantes left to assign
anno_left <- GenomicRanges::setdiff(anno_left, utr5, ignore.strand=FALSE)


cds <- GenomicRanges::intersect(anno_left,
                                anno_regions[anno_regions$type=="CDS"] ) %>%
  GenomicRanges::reduce()
anno_left <- GenomicRanges::setdiff(anno_left, cds, ignore.strand=FALSE)


intron <- GenomicRanges::intersect(anno_left,
                                   anno_regions[anno_regions$type=="intron"] ) %>%
  GenomicRanges::reduce()
anno_left <- GenomicRanges::setdiff(anno_left, intron, ignore.strand=FALSE)

# get lenght of regions for as norm factor
norm_df <- data.frame(region = c("three_prime_UTR", "five_prime_UTR", "CDS", "intron"),
                      width = c(sum(width(utr3)), sum(width(utr5)), sum(width(cds)), sum(width(intron))))
  
# normalise bs in regions by norm factor
norm_regions <- binding_sites_with_regions %>% as.data.frame(.) %>%
  group_by(region) %>%
  summarize(nts_in_BS = sum(width-2)) %>% 
  left_join(norm_df, by = "region") %>%
  mutate(norm_binding_per_reg_size = nts_in_BS / width) # mal 1000 per kb


# plot normed gene regions
ggplot(norm_regions, aes(x=region, y = norm_binding_per_reg_size))+
  geom_col()+
  theme_paper()

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

```

## Strongest binding site per gene
 
```{r}
# check strongest binding site per gene, to see if distribution between regions shifts
# --> the strongest BS per gene is more often in the 3'UTR
strongest_BS <- binding_sites_with_regions %>%
  as.data.frame(.) %>%
  group_by(gene_id) %>%
  arrange(desc(score)) %>%
  dplyr::slice(1) %>%
  ungroup() 

ggplot(strongest_BS, aes(x=region))+
  geom_bar()+
  theme_paper()

strongest_BS_norm <- strongest_BS %>%
  group_by(region)%>%
  summarize(nts_in_BS = sum(width)) %>% 
  left_join(., norm_df, by = "region") %>%
  mutate(norm_binding_per_reg_size = nts_in_BS / width)


ggplot(strongest_BS_norm, aes(x=region, y = norm_binding_per_reg_size))+
  geom_col()+
  theme_paper()

```


# Distribution of crosslinks on gene regions

```{r eval=FALSE, include=FALSE}

#overlaps of binding sites with genes
cl <- c(CL_list[[1]], CL_list[[2]])
ol_cl <- findOverlaps(anno_regions, cl, ignore.strand = FALSE, type ="any")
cl$idx <- 1:NROW(cl)

# get all possible regions overlapping with clrosslinked nucleotide
cl_2 <- cl[subjectHits(ol_cl)]
elementMetadata(cl_2) <- c(elementMetadata(cl[subjectHits(ol_cl)]), elementMetadata(anno_regions[queryHits(ol_cl)]))

# choose type highest in hiearcy
cl_3 <- as.data.frame(cl_2) %>%
  group_by(idx) %>%
  arrange(factor(type, levels = c("three_prime_UTR", "five_prime_UTR", "CDS", "intron")), .by_group = T) %>%
  dplyr::slice(1)


ggplot(cl_3, aes(x=type))+
  geom_bar()+
  stat_count(aes(label=paste0(sprintf("%1.1f", ..count../sum(..count..)*100),
                              "% \n", ..count..), y=0.5*..count..), 
             geom="text", colour="white", size=4, position=position_dodge(width=1)) +
  theme_paper()+
  ylab("Number of crosslink events")

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

```


# Meta gene profile

```{r eval=FALSE, include=FALSE}
# split for crosslinks on genes with highestbbinding in cds vs utr
BS_cds_utr <- binding_sites_with_regions[binding_sites_with_regions$region %in% c("CDS", "three_prime_UTR"),]
highest_BS_per_gene <- BS_cds_utr %>%
  as.data.frame(.) %>%
  group_by(gene_id) %>%
  arrange(desc(score), .by_group = T) %>%
  dplyr::slice(1) %>%
  ungroup() %>%
  mutate(gene_id = substr(gene_id,1,15))

# selct only crosslinks on annotated genes
idx_cl_gene_id <- findOverlaps(cl, annotation[annotation$type == "gene"])
idx_cl_gene_id <- as.data.frame(idx_cl_gene_id) %>%  filter(!duplicated(queryHits))

cl$gene_id <- NA
anno_gene <- annotation[annotation$type == "gene"]
cl[idx_cl_gene_id$queryHits]$gene_id <- anno_gene[idx_cl_gene_id$subjectHits]$gene_id

cl_on_genes <- cl[!is.na(cl$gene_id)] %>%
  as.data.frame(.) %>%
  left_join(highest_BS_per_gene[,c("gene_id", "region")], by = "gene_id") %>% 
  filter(!is.na(region))%>%
  mutate(group = region) %>%
  makeGRangesFromDataFrame(keep.extra.columns = T)

# make metaprofile 

profile <- metaGeneProfile(cl_on_genes, annotation = gff3, group="region")
profile$Plot + ylim(c(0,0.9)) + theme_paper() + theme(legend.position = NULL)

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

```

# Binding vs tpm


## RNAseq TPM vs CLIP TPM

```{r}
#########################
# Calcualte CLIP tpms
#########################

# get transcript lengths
transcript_length <- transcriptLengths(anno_txdb, with.cds_len = T)
transcript_length_mean <- transcript_length %>%
  dplyr::group_by(gene_id) %>%
  summarise(transcript_length = mean(tx_len), .groups = "keep")
transcript_length_mean$gene_id <- substr(transcript_length_mean$gene_id,1,15)
binding_sites_with_regions$gene_id <- substr(binding_sites_with_regions$gene_id,1,15)

# add mean transcript length to binding sites
binding_sites_with_regions <- binding_sites_with_regions %>% as.data.frame() %>%
  left_join(., transcript_length_mean, by = "gene_id", ) 

# add gene length for genes with no annotated transcript
binding_sites_with_regions[is.na(binding_sites_with_regions$transcript_length),]$binding_sites_with_regions <-  binding_sites_with_regions[is.na(binding_sites_with_regions$transcript_length),]$width

# tpm for all samples
binding_sites_with_regions <- binding_sites_with_regions %>% mutate(
    s1.clip.tpm = clp_rep1 / transcript_length,
    s2.clip.tpm = clp_rep2 / transcript_length,
    s3.clip.tpm = clp_rep3 / transcript_length,
    s4.clip.tpm = clp_rep4 / transcript_length
 )


# scaling factor values for samples 
scaling_factor <- binding_sites_with_regions[,c("s1.clip.tpm", "s2.clip.tpm", "s3.clip.tpm", "s4.clip.tpm")] 
scaling_factor <- colSums(scaling_factor, na.rm = T)
scaling_factor <- scaling_factor/1000000

binding_sites_with_regions <- binding_sites_with_regions %>% mutate(
    s1.clip.tpm = s1.clip.tpm / scaling_factor["s1.clip.tpm"],
    s2.clip.tpm = s2.clip.tpm / scaling_factor["s2.clip.tpm"],
    s3.clip.tpm = s3.clip.tpm / scaling_factor["s3.clip.tpm"],
    s4.clip.tpm = s4.clip.tpm  / scaling_factor["s4.clip.tpm"]) %>%
  rowwise() %>%
  mutate( mean.clip.tpm = mean(c(s1.clip.tpm, s2.clip.tpm, s3.clip.tpm, s4.clip.tpm))
    )


#############################
# add RNAseq TPM and compare
############################
binding_sites_with_regions <- binding_sites_with_regions %>% left_join(., rnaseq_counts, by = c(gene_id= "gene_id"))

ggplot(binding_sites_with_regions, aes(x = log10(mean_tpm_rnaseq), log10(mean.clip.tpm)))+
                ggrastr::rasterise(geom_point(size = 0.2), dpi = 300)+
  ggrastr::rasterise(ggpointdensity::geom_pointdensity(size = 0.2), dpi = 300)+
  geom_smooth(method = "lm", color = "black")+
  stat_cor()+
  theme_paper()
  
ggsave(filename = paste0(outpath, "tpm_comp.pdf"), width = 7, height = 5, units = "cm")


```



# Dependency of binding to tpm

```{r}
##############################
# bin binding sites over a certan RNAseq tpm cutoff
#############################

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

# make list with bs about cutoffs
RNA_tpm_cutoffs <- map(s, ~rnaseq_counts %>% filter(mean_tpm_rnaseq > .x))

# summarise number of bs over each cutoff in one data frame
RNA_tpm_cutoffs_bound <- map2_dfr(RNA_tpm_cutoffs, s, ~c(bound = nrow(.x[.x$gene_id %in% substr(binding_sites_with_regions$gene_id,1,15),]),
                                                  all_over_tpm_cut =   nrow(.x),
                                                  tpm_cut = .y))

# calculate ration of bound vs unbound for bindingsites over each cutoff
RNA_tpm_cutoffs_bound <- mutate(RNA_tpm_cutoffs_bound, ratio_bound = bound / all_over_tpm_cut)

# plot
ggplot(RNA_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), "% \n of", all_over_tpm_cut)), 
             geom="text", color= "black", size=4)+
  #ggrepel::geom_text_repel()+
  ylim(0,1)+
  theme_paper()

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

```


# Distribution of binding sites per gene

From revision:
 - the distribution of the number of sites per gene over the whole transcriptome should be shown. Is it a continuum from no binding sites to many, or is there a subset of genes intensively bound

## as histogram

```{r}
# binding_sites_with_regions <- makeGRangesFromDataFrame(binding_sites_with_regions, keep.extra.columns = T)
# bound_genes <- anno_gene[substr(anno_gene$gene_id,1,15) %in% binding_sites_with_regions$gene_id]

n_BS_per_gene <- binding_sites_with_regions %>% 
  as.data.frame() %>%
  group_by(gene_id) %>%
  summarize(n_BS = n()) 

n_BS_per_gene  <- n_BS_per_gene %>%
  mutate( n_BS_plot = case_when(n_BS > 100 ~ 100,
                                   T ~ as.numeric(n_BS)))


ggplot(n_BS_per_gene, aes(x=n_BS_plot))+
  geom_histogram(binwidth = 1)+
  theme_paper()+
  geom_vline(xintercept = mean(n_BS_per_gene$n_BS), color =  "blue")+
  geom_vline(xintercept = median(n_BS_per_gene$n_BS), color =  "green")+
  geom_vline(xintercept = quantile(n_BS_per_gene$n_BS)[2], color = "green")+
   geom_vline(xintercept = quantile(n_BS_per_gene$n_BS)[4], color =  "green")+
  xlab("Number of binding sites per gene")+
  ylab("Number of genes")

ggsave(filename = paste0(outpath, "bs_per_gene_distribution.pdf"), width = 8, height = 5, units = "cm")

```



## in bins

```{r th}
###################
# Number of BS per gene
##################

n_BS_per_gene_tpm <- left_join(n_BS_per_gene, rnaseq_counts, by = "gene_id" )


ggplot(n_BS_per_gene_tpm, aes(x=mean_tpm_rnaseq, y = n_BS))+
  geom_point()+
  #ggforce::facet_zoom(xlim=c(0,1500))+
  #geom_hline(yintercept = 25, color ="red")+
  #geom_smooth()+
  ylab("Number of binding sites per gene")+
  xlab("Gene expression [tpm]")+
  ylim(c(0,300))+
  scale_x_log10()+
  theme_paper()



n_BS_per_gene_tpm <- n_BS_per_gene_tpm %>% mutate(bin = cut(mean_tpm_rnaseq, breaks = c(-Inf,1,10,100,200,300,400,500,600,25000 )) )

median_per_bin <- n_BS_per_gene_tpm %>% group_by(bin) %>%
  summarise(median = median(n_BS))

ggplot(n_BS_per_gene_tpm, aes(x=as.factor(bin), y = n_BS))+
  # gghalves::geom_half_boxplot(outlier.size = -5, lwd = 0.5) +
  geom_boxplot()+
  #geom_point(data = median_per_bin, aes(x = bin, y = median ), color = "red")+
  theme_paper()+
  ylab("Number of binding sites per gene")+
  xlab("Binned gene expression [tpm]")+
  ylim(c(0,250))+
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

ggsave(paste0(outpath,"tpm_vs_BS_per_gene.pdf"), height = 6, width = 7, units = "cm")


```


# GO

```{r eval=FALSE, include=FALSE}
# background: all genes with crosslinks
expr_genes <- subsetByOverlaps(annotation_gene, c(CL_list[[1]],CL_list[[2]]))

# msigdb_available(species = "Homo sapiens")
reactome_geneset <- msigdb_gsets("Homo sapiens", "C2", "CP:REACTOME")
cell_comp_geneset <- msigdb_gsets("Homo sapiens", "C5", "CC")


# background
bg <- unique(expr_genes$gene_name)

# Test USS groups vs gene sets
hyps_reactome <-  hypeR(bound_genes_gg$gene_name, reactome_geneset, test="hypergeometric", fdr=1, pval = 1, background = bg )
hyps_cell_comp <- hypeR(bound_genes_gg$gene_name, cell_comp_geneset, test="hypergeometric", fdr=1, pval = 1, background = bg) 



# plot reactome
reactome_df <- hyps_reactome$data %>% mutate(label = tolower(label), generatio = overlap / geneset)

ggplot(reactome_df[1:25,], aes(y = generatio, x = factor(label, level = rev(unique(label)))))+
  geom_col()+
  coord_flip()+
  theme_paper()+
  scale_size(range = c(1, 3))

xlsx::write.xlsx(reactome_df, paste0(outpath, "iCLIP_reactome_list.xlsx"))
ggsave(paste0(outpath, Sys.Date(), "reactome_iCLIP_paper.pdf"), width = 17, height = 8, units = "cm")


# plot cellular components 
cc_df <- hyps_cell_comp$data %>% mutate(label = tolower(label), generatio = overlap / geneset)

ggplot(cc_df[1:25,], aes(y = generatio, x = factor(label, level = rev(unique(label)))))+
  geom_col()+
  coord_flip()+
  theme_paper()+
  scale_size(range = c(1, 3))

xlsx::write.xlsx(cc_df, paste0(outpath, "iCLIP_GO_cc_list.xlsx"))
ggsave(paste0(outpath, Sys.Date(), "cc_iCLIP_paper.pdf"), width = 17, height = 8, units = "cm")


# cellular components connected to neurons
dendritic_cc <- cc_df[grepl(cc_df$label, pattern="dendritic"),]
neuron_cc <- cc_df[grepl(cc_df$label, pattern="neuron"),]
axon_cc <- cc_df[grepl(cc_df$label, pattern="axon"),]

neuro_dendritic_cc <- rbind(dendritic_cc, neuron_cc, axon_cc)

ggplot(neuro_dendritic_cc[neuro_dendritic_cc$fdr < 0.05,], aes(y = generatio, x = factor(label, level = rev(unique(label))), fill = -log10(fdr) ))+
  geom_col()+
  coord_flip()+
  theme_paper()+
  scale_size(range = c(1, 3))

ggsave(paste0(outpath, Sys.Date(), "cc_dndritic_iCLIP_paper.pdf"), width = 17, height = 8, units = "cm")


```


# Output
- binding sites with additional infos (binding_sites_characterized.rds)

```{r eval=FALSE, include=FALSE}

saveRDS(binding_sites_with_regions,  "/Users/melinaklostermann/Documents/projects/PURA/Molitor-et-al-2022/binding_sites_characterized.rds")

```