---
title: "Figure_2_code"
author: "Alex Nesta"
date: "2023-09-22"
output: html_document
---

## Load Required Packages and set working paths

```{r}
## Packages
library(tidyverse)
library(rtracklayer)
library(plyranges)
library(scales)
library(ggpubr)
library(ggVennDiagram)

## Code Path
code_path <- ("/Volumes/Untitled/code_for_figures/")
setwd(code_path)

## Data Path
data_path <- ("/Volumes/Untitled/data_for_figures/")

## Plots Path
plots_path <- ("/Volumes/Untitled/plots_for_figures/")

## Tables Path
tables_path <- ("/Volumes/Untitled/tables_for_figures/")

```

## Load Data

```{r}

source("./extract_suppa_junctions_v2.R")

repeats <- read_tsv(paste0(data_path, "rmsk_simp.txt"))

# Working with suppa data now
new_df <- read_tsv(paste0(data_path, "all_clusters_differentialSplicing.tsv"))

```

# Process repeats
```{r}


rmsk <- repeats %>% dplyr::select(chr = genoName, start = genoStart, end = genoEnd, strand, repName, repClass, repFamily)
rmsk.gr <- as(rmsk, "GRanges")

#We need to preserce chr, start, and end for downstream filtering
simp_rmsk <- rmsk.gr %>% 
  filter(repClass %in% c("LTR", "LINE", "SINE", "DNA")) %>%
  mutate(repChr = seqnames(.)) %>%
  mutate(repStart = start(.)) %>%
  mutate(repEnd = end(.)) %>%
  mutate(repStrand = strand)

export.bed(simp_rmsk, "simp_rmsk.bed")


```


# Process suppa data
```{r}
colnames(new_df)

#Remove _250bp from type
new_df <- new_df %>% separate(type, sep = "_", into = "type")


# Extracted Splice Events
extracted <- extract_suppa(new_df)

colnames(extracted)

## Backup critical columns before the intersect
grange_prep <- extracted %>% 
  mutate(Junction = junction_bak, start = start_1, end = end_1, start_1_bak = start_1, end_1_bak = end_1, start_2_bak = start_2, end_2_bak = end_2)

#Save the modified DF for other analyses
grange_prep %>% write_tsv(paste0(tables_path, "as_events_coords_extracted.tsv"))

```

## Find TE overlapping splicing events

```{r}


# Make a grange object so we can do intersects
diff_exons <- makeGRangesFromDataFrame(grange_prep,
                         seqnames.field = "chr",
                         start.field = "start",
                         end.field = "end",
                         strand.field = "strand",
                         keep.extra.columns = TRUE)

diff_exons <- diff_exons %>% mutate(chr_orig = seqnames(diff_exons), 
                                    start_orig = start(diff_exons), 
                                    end_orig = end(diff_exons))

### Just get the 5' or 3' ends of each exon

diff_exons_5p <- mutate(anchor_5p(diff_exons), width = 1)
diff_exons_3p <- mutate(anchor_3p(diff_exons), width = 1)

#Intersect and label 5' and 3' ends of exons with repeat masker
overlaps_5p <- find_overlaps(diff_exons_5p, simp_rmsk)

overlaps_3p <- find_overlaps(diff_exons_3p, simp_rmsk)

overlaps_5p_df <- as_tibble(overlaps_5p %>% mutate(TE_5p = TRUE)) 
overlaps_3p_df <- as_tibble(overlaps_3p %>% mutate(TE_3p = TRUE))

all_hits <- bind_rows(overlaps_5p_df, overlaps_3p_df) 

#MAIN FILTERING CRITERIA. NOTE: (& != &&)

passing_events <- all_hits %>% mutate(keep = case_when(event_type == "AF" & TE_5p == TRUE ~ TRUE,
                                                       event_type == "AL" & TE_3p == TRUE ~ TRUE,
                                                       event_type == "AL" & TE_5p == TRUE ~ TRUE,
                                                       event_type == "SE" & TE_5p == TRUE ~ TRUE,
                                                       event_type == "A5" & TE_3p == TRUE ~ TRUE,
                                                       event_type == "A3" & TE_5p == TRUE ~ TRUE,
                                                       event_type == "MX" & TE_5p == TRUE ~ TRUE))

passing_events_ref_final <-  
  passing_events %>%
  filter(keep == TRUE) %>%
  mutate(TE_IN_REF = TRUE)



########## SPECIAL OVERLAP FOR ALTERNATIVE FIRST AND ALTERNATIVE LAST SECONDARY TEs.
# Make a grange object so we can do intersects
diff_exons_alt <- grange_prep %>%
  filter(event_type %in% c("AF", "AL")) %>%
  dplyr::select(-start, -end) %>%
  makeGRangesFromDataFrame(seqnames.field = "chr",
                           start.field = "start_2_bak",
                           end.field = "end_2_bak",
                           strand.field = "strand",
                           keep.extra.columns = TRUE)

diff_exons_alt <- diff_exons_alt %>% mutate(chr_orig = seqnames(diff_exons_alt), 
                                            start_2_bak = start(diff_exons_alt), 
                                            end_2_bak = end(diff_exons_alt))

### Just get the 5' or 3' ends of each exon

diff_exons_5p_alt <- mutate(anchor_5p(diff_exons_alt), width = 1)
diff_exons_3p_alt <- mutate(anchor_3p(diff_exons_alt), width = 1)

#Intersect and label 5' and 3' ends of exons with repeat masker
overlaps_5p_alt <- find_overlaps(diff_exons_5p_alt, simp_rmsk)
overlaps_3p_alt <- find_overlaps(diff_exons_3p_alt, simp_rmsk)

overlaps_5p_df_alt <- as_tibble(overlaps_5p_alt %>% mutate(TE_5p = TRUE)) 
overlaps_3p_df_alt <- as_tibble(overlaps_3p_alt %>% mutate(TE_3p = TRUE))

all_hits_alt <- bind_rows(overlaps_5p_df_alt, overlaps_3p_df_alt) 

passing_events_alt <- all_hits_alt %>% mutate(keep = case_when(event_type == "AF" & TE_5p == TRUE ~ TRUE,
                                                           event_type == "AL" & TE_3p == TRUE ~ TRUE,
                                                           event_type == "AL" & TE_5p == TRUE ~ TRUE))

passing_events_alt_final <- passing_events_alt %>%
  filter(keep == TRUE) %>%
  mutate(TE_IN_ALT = TRUE)

######### Special intersect just for retained introns
retained_introns <- diff_exons %>% filter(event_type == "RI")

RI_overlaps <- find_overlaps(retained_introns, simp_rmsk) %>% mutate(keep = TRUE) %>% as_tibble()

#because retained_introns overlap multiple repeats, the repeats need to be collapsed.

RI_overlaps_collapsed <- RI_overlaps %>% 
  group_by(across(-starts_with("rep"))) %>%
  summarise(repName = paste(repName, collapse =","),
            repClass = paste(repClass, collapse = ","),
            repFamily = paste(repFamily, collapse = ","),
            repChr = paste(repChr, collapse = ","),
            repStart = paste(repStart, collapse = ","),
            repEnd = paste(repEnd, collapse = ","),
            repStrand = paste(repStrand, collapse = ",")) %>% 
  ungroup() %>%
  mutate(keep = TRUE)


passing_RI_events_ref_final <-  
  RI_overlaps_collapsed %>%
  filter(keep == TRUE) %>%
  mutate(TE_IN_REF = TRUE)

            


RI_char <- as_tibble(lapply(passing_RI_events_ref_final, as.character))
pass_char <- as_tibble(lapply(passing_events_ref_final, as.character))
pass_alt_char <- as_tibble(lapply(passing_events_alt_final, as.character)) %>% distinct(junction_bak, .keep_all = T)

RI_and_pass_ref <- bind_rows(pass_char, RI_char) %>% distinct(junction_bak, .keep_all = T)

RI_and_pass_ref %>% write_tsv(paste0(tables_path, "TE_in_ref.tsv"))
pass_alt_char %>% write_tsv(paste0(tables_path, "TE_in_alt.tsv"))

higher_te_cancer_ref <- RI_and_pass_ref %>% filter(deltaPSI_TvsN.x > 0)
higher_te_normal_ref <- RI_and_pass_ref %>% filter(deltaPSI_TvsN.x < 0)

higher_te_cancer_alt <- pass_alt_char %>% filter(deltaPSI_TvsN.x < 0)  
higher_te_normal_alt <- pass_alt_char %>% filter(deltaPSI_TvsN.x > 0) 


##### TOTAL TES IN CANCER AND TOTAL TES IN NORMAL
# total_te_cancer <- bind_rows(higher_te_cancer_ref, higher_te_cancer_alt)
# total_te_normal <- bind_rows(higher_te_normal_ref, higher_te_normal_alt)
total_te_cancer <- higher_te_cancer_ref
total_te_normal <- higher_te_normal_ref

total_te_cancer %>% write_tsv(paste0(tables_path, "TE_higher_cancer.tsv"))
total_te_normal %>% write_tsv(paste0(tables_path, "TE_higher_normal.tsv"))
all_events_combined <- bind_rows(total_te_cancer, total_te_normal)
all_events_combined %>% write_tsv(paste0(tables_path, "all_TE_events.tsv"))

```





# Do a pathway enrichment for the TE higher in cancer or te higher in normal events

```{r}


library(enrichR)

alldbs <- listEnrichrDbs()

dbs <- c("KEGG_2019_Human", "Reactome_2016", "GO_Biological_Process_2018", "GO_Cellular_Component_2018", "GO_Molecular_Function_2018", "MSigDB_Hallmark_2020", "MSigDB_Oncogenic_Signatures")

# Get the gene names for the TE higher in cancer and TE higher in normal events
te_cancer_genes <- total_te_cancer %>% pull(Gene) %>% unique()
te_normal_genes <- total_te_normal %>% pull(Gene) %>% unique()

enriched_cancer <- enrichr(te_cancer_genes, dbs)
enriched_normal <- enrichr(te_normal_genes, dbs)

plotEnrich(enriched_cancer[[7]], showTerms = 20, numChar = 40, y = "Count", orderBy = "P.value")
plotEnrich(enriched_normal[[7]], showTerms = 20, numChar = 40, y = "Count", orderBy = "P.value")


plotEnrich(enriched_cancer[[6]], showTerms = 20, numChar = 40, y = "Count", orderBy = "P.value")
plotEnrich(enriched_normal[[6]], showTerms = 20, numChar = 40, y = "Count", orderBy = "P.value")


# Make a venn diagram of the genes
library(ggvenn)
venn_genes <- list(cancer_genes = te_cancer_genes, normal_genes = te_normal_genes)
venn_genes_plot <- ggvenn(venn_genes)


library(patchwork)
# Save plot of venn diagram genes with an AS TE higher in normal or cancer
pdf(paste0(plots_path, "venn_genes_TE_cancer_normal.pdf"))
venn_genes_plot
dev.off()

```




# Check for breast cancer subtype enriched events

```{r}

all_events_combined


colnames(all_events_combined)
#  [61] "Basal_Like_log10_pval_adj"  "LuminalA_log10_pval_adj"    "LuminalB_log10_pval_adj"  "HER2_log10_pval_adj"   
# For the TE higher cancer events, check for breast cancer subtype enrichment

# First select those higher in cancer
subtype_enrichment <-
all_events_combined %>%
  filter(event_type %in% c("AF", "AL", "SE", "RI")) %>%
    mutate(pos_neg = case_when(TE_IN_REF == "TRUE" & deltaPSI_TvsN.x > 0 ~ "higher in tumor",
                               # TE_IN_ALT == "TRUE" & deltaPSI_TvsN.x < 0 ~ "higher in tumor",
                               TE_IN_REF == "TRUE" & deltaPSI_TvsN.x < 0 ~ "higher in normal",
                               # TE_IN_ALT == "TRUE" & deltaPSI_TvsN.x > 0 ~ "higher in normal"
                               )) %>%
    filter(pos_neg == "higher in tumor") %>%
    filter(Basal_Like_log10_pval_adj > 1.3 | LuminalA_log10_pval_adj > 1.3 | LuminalB_log10_pval_adj > 1.3 | HER2_log10_pval_adj > 1.3) %>%

    # select(junction_bak, Gene, event_type, TE_IN_REF, TE_IN_ALT, Counts.x, absDeltaPSI_TvsN.x, Basal_Like_log10_pval_adj, LuminalA_log10_pval_adj, LuminalB_log10_pval_adj, HER2_log10_pval_adj) %>%
    write_tsv(paste0(tables_path, "TE_higher_cancer_subtype_enrichment.tsv"))


subtype_enrichment %>% distinct(Gene)

# Show group by breast cancer subtype and show top 5 most common events (by Counts.x) for each subtype and fill by absdeltaPSI_TvsN.x
top_events_by_subtype <-
  subtype_enrichment %>%
  mutate_at(c("Basal_Like_log10_pval_adj", "LuminalA_log10_pval_adj", "LuminalB_log10_pval_adj", "HER2_log10_pval_adj", "absDeltaPSI_TvsN.x", "Counts.x"), as.numeric) %>%
  pivot_longer(cols = c(Basal_Like_log10_pval_adj, LuminalA_log10_pval_adj, LuminalB_log10_pval_adj, HER2_log10_pval_adj), names_to = "subtype", values_to = "subtype_log10_pval_adj") %>%
  filter(subtype_log10_pval_adj > 1.3) %>%
  # Group by subtype and gene and event type. Select only the top pval for each
  group_by(subtype, Gene, event_type) %>%
  dplyr::slice(1) %>%
  group_by(subtype) %>%
  dplyr::slice(1:8) %>%
  mutate(top_subtype_for_event = subtype,
         top_pval_for_event = subtype_log10_pval_adj) %>%
  mutate(gene_type = paste0(Gene,"-",event_type),
         gene_type = forcats::fct_reorder(gene_type, subtype_log10_pval_adj)) %>%
  ungroup()
           

top_events_by_subtype %>% 
  select(Gene, event_type, Survival_overall_sig)


top_event_annotation <- top_events_by_subtype %>% 
  select(junction_bak, top_subtype_for_event, top_pval_for_event)


library(patchwork)  

# Plot top 5 events by subtype and facet by subtype
pdf(paste0(plots_path, "TE-mediated-splice-class-diverged-subtype-enrichment.pdf"), height = 7, width = 4)
p1 <-
  top_events_by_subtype %>%
  mutate(event_type = if_else(event_type == "CE", "CA", event_type)) %>%
      mutate(subtype = case_when(subtype == "HER2_log10_pval_adj" ~ "HER2",
                                           subtype == "Basal_Like_log10_pval_adj" ~ "Basal",
                                           subtype == "LuminalB_log10_pval_adj" ~ "Luminal B",
                                           subtype == "LuminalA_log10_pval_adj" ~ "Luminal A")) %>%
  ggplot(aes(x = gene_type, y = Counts.x, fill = absDeltaPSI_TvsN.x)) +
  geom_col() +
  facet_wrap(~subtype, ncol = 1, scales = "free_y") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5)) +
  labs(x = "Gene Symbol - Splice event", y = "Number of tumors", fill = "Percent Spliced In") +
  # theme(legend.position = "right") + 
  coord_flip() +
  geom_hline(aes(yintercept = 50), 
             color = "red") +
  scale_fill_gradient2(limits = c(-1.0, 1.0), breaks=c(1.0, 0.5, 0, -0.5, -1.0), low = "#00478F", mid="white", high = "#FF5D00", midpoint = 0, na.value = NA)


preproc <- 
subtype_enrichment %>%
  filter(junction_bak %in% top_events_by_subtype$junction_bak) %>%
  mutate(gene_type = paste0(Gene,"-",event_type)) %>%
    mutate_at(c("Basal_Like_log10_pval_adj", "LuminalA_log10_pval_adj", "LuminalB_log10_pval_adj", "HER2_log10_pval_adj", "absDeltaPSI_TvsN.x", "Counts.x"), as.numeric) %>%
  pivot_longer(cols = c(Basal_Like_log10_pval_adj, LuminalA_log10_pval_adj, LuminalB_log10_pval_adj, HER2_log10_pval_adj), names_to = "subtype", values_to = "max_subtype_log10_pval_adj") %>%
    mutate(subtype = case_when(
    subtype == "Basal_Like_log10_pval_adj" ~ "Basal", 
    subtype == "HER2_log10_pval_adj" ~ "HER2", 
    subtype == "LuminalA_log10_pval_adj" ~ "LumA", 
    subtype == "LuminalB_log10_pval_adj" ~ "LumB")
    ) %>%
  select(junction_bak, gene_type, subtype, max_subtype_log10_pval_adj)

p2 <-
  preproc %>%
  left_join(subtype_enrichment) %>%
  left_join(top_event_annotation) %>%
  mutate(top_subtype_for_event = case_when(top_subtype_for_event == "HER2_log10_pval_adj" ~ "HER2",
                                           top_subtype_for_event == "Basal_Like_log10_pval_adj" ~ "Basal",
                                           top_subtype_for_event == "LuminalB_log10_pval_adj" ~ "Luminal B",
                                           top_subtype_for_event == "LuminalA_log10_pval_adj" ~ "Luminal A")) %>%
  group_by(top_subtype_for_event) %>%
  mutate(gene_type = forcats::fct_reorder(gene_type, top_pval_for_event)) %>%
  ggplot(aes(x = gene_type, y = subtype, fill = max_subtype_log10_pval_adj)) +
  geom_tile() +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5)) +
  labs(x = "Gene Symbol - Splice event", y = "Breast cancer subtype", fill = "-log10(pVal)") +
  # theme(legend.position = "right") + 
  coord_flip() +
  facet_wrap(~top_subtype_for_event, ncol = 1, scales = "free_y")+
    theme(
        axis.text.y = element_blank(),
        axis.title.y = element_blank(),
        axis.ticks.y = element_blank())

p1 + p2 +
  plot_layout(guides = "collect",
                          # axes = "collect",
                          # tag_level = "new",
                          # axis_titles = "collect",
              ) & theme(              legend.position = "bottom", legend.box="vertical"
) #&guides(fill=guide_legend(nrow=2,byrow=TRUE))

dev.off()




```




# Survival events
```{r}

library(ggridges)

te_surv_events <-
  all_events_combined %>%
    mutate(pos_neg = case_when(TE_IN_REF == "TRUE" & deltaPSI_TvsN.x > 0 ~ "higher in tumor",
                               # TE_IN_ALT == "TRUE" & deltaPSI_TvsN.x < 0 ~ "higher in tumor",
                               TE_IN_REF == "TRUE" & deltaPSI_TvsN.x < 0 ~ "higher in normal",
                               # TE_IN_ALT == "TRUE" & deltaPSI_TvsN.x > 0 ~ "higher in normal"
                               )) %>%
    filter(pos_neg == "higher in tumor") %>%
  filter(event_type %in% c("AF", "CA" = "SE", "AL", "RI")) %>%
  filter(Survival_global_pval < 0.05)

te_surv_events %>%
  openxlsx::write.xlsx(file = paste0(plots_path, "TE-mediated-splice-class-diverged-survival.xlsx"), as.data.frame = TRUE)



significant_survival_events <-
all_events_combined %>%
    mutate(pos_neg = case_when(TE_IN_REF == "TRUE" & deltaPSI_TvsN.x > 0 ~ "higher in tumor",
                               # TE_IN_ALT == "TRUE" & deltaPSI_TvsN.x < 0 ~ "higher in tumor",
                               TE_IN_REF == "TRUE" & deltaPSI_TvsN.x < 0 ~ "higher in normal",
                               # TE_IN_ALT == "TRUE" & deltaPSI_TvsN.x > 0 ~ "higher in normal"
                               )) %>%
    filter(pos_neg == "higher in tumor") %>%
  filter(event_type %in% c("AF", "SE", "AL", "RI")) %>%
  filter(Survival_global_pval < 0.05)


# plot the set of splicing events that are significant in the survival analysis

# Make a bar plot with number of patients with event
p1 <- significant_survival_events %>%
  ungroup() %>%
  mutate(Counts.x = as.numeric(Counts.x)) %>%
  mutate(Gene = forcats::fct_reorder(Gene, Counts.x)) %>%
  mutate(event_type = factor(event_type, levels = c("AF", "CA" = "SE", "AL", "RI")))%>%
  ggplot(aes(x = Gene, y = Counts.x, fill = -log10(as.numeric(Survival_global_pval)))) +
  geom_col() +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5)) +
  labs(x = "Gene Symbol", y = "Number of patients with event", fill = "-log10(p) Survival") +
  coord_flip() +
  facet_wrap(~event_type, ncol = 1, scales = "free_y") +
  ggtitle("Survival associated events")

p1
# Make a heatmap with ggplot showing repClass
p2 <- significant_survival_events %>%
  ungroup() %>%
  mutate(Counts.x = as.numeric(Counts.x)) %>%
  mutate(Gene = forcats::fct_reorder(Gene, Counts.x)) %>%
  mutate(event_type = factor(event_type, levels = c("AF","CA" = "SE", "AL", "RI")))%>%
  ggplot(aes(x = Gene, y = event_type, fill = repFamily)) +
  geom_tile(color="black") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5)) +
  labs(x = "Gene Symbol", y = "Event type", fill = "Repeat Family") +
  coord_flip() +
  facet_wrap(~event_type, ncol = 1, scales = "free") +
  scale_fill_manual(values = c("Alu" = "darkorange", "SVA" = "darkgreen", "L1" = "darkblue", "L2" = "lightblue", "ERVL-MaLR" = "darkred", "MIR" = "lightgreen", "RTE-X" = "forestgreen", "MER117,MER21B,MIR,MIR" = "black")) +
  ggtitle("TE\nFamily")+
  # Hide the x axis and y axis labels
  theme(axis.text.x = element_blank(),
        axis.title.x = element_blank(),
        axis.text.y = element_blank(),
        axis.title.y = element_blank(),
        axis.ticks = element_blank())

p2
# Do it for repeat class now
p3 <- significant_survival_events %>%
  ungroup() %>%
  mutate(Counts.x = as.numeric(Counts.x)) %>%
  mutate(Gene = forcats::fct_reorder(Gene, Counts.x)) %>%
  mutate(event_type = factor(event_type, levels = c("AF", "CA" = "SE", "AL", "RI")))%>%
  filter(event_type == "RI") %>%
  distinct(Gene, repClass, event_type, repFamily, repName, repFamily, repClass, repStrand) %>%
  ggplot(aes(x = Gene, y = event_type, fill = repClass)) +
  geom_tile(color="black") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5)) +
  labs(x = "Gene Symbol", y = "Event type", fill = "Repeat class") +
  coord_flip() +
  facet_wrap(~event_type, ncol = 1, scales = "free") +
  scale_fill_manual(values = c("SINE" = "orange", "LINE" = "blue", "LTR" = "red", "DNA" = "green", "Unknown" = "gray")) +
  ggtitle("TE\nClass")+
  # Hide the x axis and y axis labels
  theme(axis.text.x = element_blank(),
        axis.title.x = element_blank(),
        axis.text.y = element_blank(),
        axis.title.y = element_blank(),
        axis.ticks = element_blank())
p3

# Repeat orientation
p4 <- significant_survival_events %>%
  ungroup() %>%
    mutate("repOrientation" = (case_when(
                                    strand == repStrand ~ "sense",
                                    strand != repStrand ~ "antisense", 
                                    T ~ "Multiple"))) %>%
  mutate(Counts.x = as.numeric(Counts.x)) %>%
  mutate(Gene = forcats::fct_reorder(Gene, Counts.x)) %>%
  mutate(event_type = factor(event_type, levels = c("AF","CA" = "SE", "AL", "RI")))%>%
  ggplot(aes(x = Gene, y = event_type, fill = repOrientation)) +
  geom_tile(color="black") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5)) +
  labs(x = "Gene Symbol", y = "Event type", fill = "Repeat orientation") +
  coord_flip() +
  facet_wrap(~event_type, ncol = 1, scales = "free") +
  scale_fill_manual(values = c("sense" = "#D8E1E7", "antisense" = "#606060", "Multiple" = "black", "unknown" = "gray")) +
  ggtitle("TE\norientation") +
  # Hide the x axis and y axis labels
  theme(axis.text.x = element_blank(),
        axis.title.x = element_blank(),
        axis.text.y = element_blank(),
        axis.title.y = element_blank(),
        axis.ticks = element_blank())
p4

# Plot the PSI values
p5 <- significant_survival_events %>%
  ungroup() %>%
  mutate(Counts.x = as.numeric(Counts.x)) %>%
  mutate(Gene = forcats::fct_reorder(Gene, Counts.x)) %>%
  mutate(event_type = factor(event_type, levels = c("AF","CA" = "SE", "AL", "RI")))%>%
  # Pivot longer everything that starts with mean_PSI
  # Make all mean_PSI columns numeric
  mutate_at(vars(meanPSI_T.x, starts_with("mean_PSI")), as.numeric) %>%
  # Select only the columns that start with mean_PSI
  select(Gene, event_type, meanPSI_T.x, starts_with("mean_PSI")) %>%
  select(-mean_PSI_controls, -mean_PSI_gtex_all) %>%
  pivot_longer(cols = c(meanPSI_T.x, starts_with("mean_PSI")), names_to = "Tissue", values_to = "PSI_value") %>%
  # Conditionally rename the Tissue column
  mutate(Tissue = case_when(
    Tissue == "meanPSI_T.x" ~ "Breast tumor TCGA",
    Tissue == "mean_PSI_adjacent_TCGA" ~ "Breast adjacent TCGA",
    Tissue == "mean_PSI_gtex_breast" ~ "GTEx breast",
    Tissue == "mean_PSI_gtex_adipose" ~ "GTEx adipose",
    Tissue == "mean_PSI_gtex_blood" ~ "GTEx blood",
    Tissue == "mean_PSI_gtex_blood_vessel" ~ "GTEx blood vessel",
    Tissue == "mean_PSI_gtex_brain" ~ "GTEx brain",
    Tissue == "mean_PSI_gtex_colon" ~ "GTEx colon",
    Tissue == "mean_PSI_gtex_esophagus" ~ "GTEx esophagus",
    Tissue == "mean_PSI_gtex_heart" ~ "GTEx heart",
    Tissue == "mean_PSI_gtex_kidney" ~ "GTEx kidney",
    Tissue == "mean_PSI_gtex_liver" ~ "GTEx liver",
    Tissue == "mean_PSI_gtex_lung" ~ "GTEx lung",
    Tissue == "mean_PSI_gtex_muscle" ~ "GTEx muscle",
    Tissue == "mean_PSI_gtex_nerve" ~ "GTEx nerve",
    Tissue == "mean_PSI_gtex_ovary" ~ "GTEx ovary",
    Tissue == "mean_PSI_gtex_pancreas" ~ "GTEx pancreas",
    Tissue == "mean_PSI_gtex_pituitary" ~ "GTEx pituitary",
    Tissue == "mean_PSI_gtex_prostate" ~ "GTEx prostate",
    Tissue == "mean_PSI_gtex_skin" ~ "GTEx skin",
    Tissue == "mean_PSI_gtex_small_intestine" ~ "GTEx small intestine",
    Tissue == "mean_PSI_gtex_spleen" ~ "GTEx spleen",
    Tissue == "mean_PSI_gtex_stomach" ~ "GTEx stomach",
    Tissue == "mean_PSI_gtex_testis" ~ "GTEx testis",
    Tissue == "mean_PSI_gtex_thyroid" ~ "GTEx thyroid",
    Tissue == "mean_PSI_gtex_uterus" ~ "GTEx uterus",
    TRUE ~ Tissue
  )) %>%
  # Make sure the first three columns are Breast tumor TCGA, Breast adjacent TCGA, GTEx breast, the rest are alphabetically ordered
mutate(Tissue = factor(Tissue, levels = unique(c("Breast tumor TCGA", "Breast adjacent TCGA", "GTEx breast", sort(setdiff(unique(Tissue), c("Breast tumor TCGA", "Breast adjacent TCGA", "GTEx breast"))))))) %>%
  ggplot(aes(x = Gene, y = Tissue, fill = PSI_value)) +
  geom_tile(color="black") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5)) +
  labs(x = "Gene Symbol", y = "Event type", fill = "PSI value") +
  coord_flip() +
  facet_wrap(~event_type, ncol = 1, scales = "free_y") +
  scale_fill_viridis_c() +
  ggtitle("PSI values across tissues") +
  # Hide the x axis and y axis labels
  theme(axis.text.y = element_blank(),
        axis.title.y = element_blank(),
        axis.ticks.y = element_blank(),
        axis.title.x = element_blank()) 
# p5

# Make a plot for subtype enrichment
p9 <- 
significant_survival_events %>%
  ungroup() %>%
  mutate(Counts.x = as.numeric(Counts.x)) %>%
  mutate(Gene = forcats::fct_reorder(Gene, Counts.x)) %>%
  mutate(event_type = factor(event_type, levels = c("AF","CA" = "SE", "AL", "RI")))%>%
  select(Gene, event_type, ends_with("_log10_pval_adj")) %>%
  # pivot longer everything that starts with subtype_enrichment
  pivot_longer(cols = ends_with("_log10_pval_adj"), names_to = "subtype", values_to = "subtype_enrichment") %>%
  mutate(subtype_enrichment = as.numeric(subtype_enrichment),
         # Rename subtypes to plot easier
         subtype = case_when( subtype == "Basal_Like_log10_pval_adj" ~ "Basal",
                              subtype == "HER2_log10_pval_adj" ~ "HER2",
                              subtype == "LuminalA_log10_pval_adj" ~ "Luminal A",
                              subtype == "LuminalB_log10_pval_adj" ~ "Luminal B")
         # subtype_enrichment = if_else(subtype_enrichment < 1.3, NA, subtype_enrichment)
         ) %>%
  ggplot(aes(x = Gene, y = subtype, fill = subtype_enrichment)) +
  geom_tile(color="black") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5)) +
  labs(x = "Gene Symbol", y = "Event type", fill = "Subtype enrichment\n-log10(pval) enrichment") +
  coord_flip() +
  facet_wrap(~event_type, ncol = 1, scales = "free_y") +
  ggtitle("Subtype\nenrichment") +
  # Change the fill scale from gray to white to blue and make the center color 1.3
  scale_fill_gradient2(low = "white", mid = "gray", high = "firebrick", midpoint = 1.3) +
  # Hide the x axis and y axis labels
  theme(axis.text.y = element_blank(),
        axis.title.y = element_blank(),
        axis.ticks.y = element_blank(),
        axis.title.x = element_blank())


# Make a plot for transcriptome origin
p6 <-  significant_survival_events %>%
  ungroup() %>%
  mutate(Counts.x = as.numeric(Counts.x)) %>%
  mutate(Gene = forcats::fct_reorder(Gene, Counts.x)) %>%
  mutate(event_type = factor(event_type, levels = c("AF", "CA" = "SE", "AL", "RI")))%>%
  ggplot(aes(x = Gene, y = event_type, fill = transcriptome_origin)) +
  geom_tile(color="black") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5)) +
  labs(x = "Gene Symbol", y = "Event type", fill = "Isoform detection") +
  coord_flip() +
  facet_wrap(~event_type, ncol = 1, scales = "free") +
  scale_fill_manual(values = c("GENCODE" = "gray", "LR-seq" = "purple", "unknown" = "gray")) +
  ggtitle("Isoform\ndetection") +
  # Hide the x axis and y axis labels
  theme(axis.text.x = element_blank(),
        axis.title.x = element_blank(),
        axis.text.y = element_blank(),
        axis.title.y = element_blank(),
        axis.ticks = element_blank())

# p6

# Make a plot for survival prognosis
p7 <- significant_survival_events %>%
  ungroup() %>%
  mutate(Counts.x = as.numeric(Counts.x)) %>%
  mutate(Gene = forcats::fct_reorder(Gene, Counts.x)) %>%
  mutate(event_type = factor(event_type, levels = c("AF","CA" = "SE", "AL", "RI")))%>%
  ggplot(aes(x = Gene, y = event_type, fill = Survival_prognosis)) +
  geom_tile(color="black") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5)) +
  labs(x = "Gene Symbol", y = "Event type", fill = "Survival prognosis") +
  coord_flip() +
  facet_wrap(~event_type, ncol = 1, scales = "free") +
  scale_fill_manual(values = c("Favorable" = "green", "Unfavorable" = "red")) +
  ggtitle("Survival\nprognosis") +
  # Hide the x axis and y axis labels
  theme(axis.text.x = element_blank(),
        axis.title.x = element_blank(),
        axis.text.y = element_blank(),
        axis.title.y = element_blank(),
        axis.ticks = element_blank())
# p7

# plot TE_IN_REF or TE_IN_ALT
p8 <- significant_survival_events %>%
  ungroup() %>%
  mutate(Counts.x = as.numeric(Counts.x)) %>%
  mutate(Gene = forcats::fct_reorder(Gene, Counts.x)) %>%
  mutate(event_type = factor(event_type, levels = c("AF", "CA" = "SE", "AL", "RI")))%>%
  # select relevant columns and pivot long for ggplot
  select(Gene, event_type, TE_IN_REF) %>%
  mutate(TE_IN_REF = if_else(TE_IN_REF == T, "TE +PSI", "TE -PSI", missing = "TE -PSI")) %>%
  ggplot(aes(x = Gene, y = event_type, fill = TE_IN_REF)) +
  geom_tile(color="black") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5)) +
  labs(x = "Gene Symbol", y = "Event type", fill = "TE in reference") +
  coord_flip() +
  facet_wrap(~event_type, ncol = 1, scales = "free") +
  scale_fill_manual(values = c("TE +PSI" = "yellow", "TE -PSI" = "blue")) +
  ggtitle("TE +/-\nREF") +
  # Hide the x axis and y axis labels
  theme(axis.text.x = element_blank(),
        axis.title.x = element_blank(),
        axis.text.y = element_blank(),
        axis.title.y = element_blank(),
        axis.ticks = element_blank())
pdf(paste0(plots_path, "TE-mediated-splice-class-diverged-survival_no_legend.pdf"), height = 7.25, width = 10)
  p1+p7+p2+p4+p5+p9+p6+
  plot_layout(guides = "collect",
              axes = "collect",
              axis_titles = "collect",
              widths = c(1.5, .5, .5, .5, 2,1,.5,.5)
  ) & 
    theme(legend.position = "null", 
          legend.box = "horizontal",
          legend.margin = margin(t = 5, b = 5),
          legend.spacing.x = unit(1, "cm"),
          legend.key.size = unit(0.8, "lines"),
          legend.text = element_text(size = 7)) &
    guides(
      fill = guide_legend(ncol = 2, byrow = FALSE),
      color = guide_legend(ncol = 2, byrow = FALSE)
    )
dev.off()



pdf(paste0(plots_path, "TE-mediated-splice-class-diverged-survival_legend.pdf"), height = 10, width = 12)
  p1+p7+p2+p4+p5+p9+p6+
  plot_layout(guides = "collect",
              axes = "collect",
              axis_titles = "collect",
              widths = c(1.5, .5, .5, .5, 2,1,.5,.5)
  ) & 
    theme(legend.position = "bottom", 
          legend.box = "vertical",
          legend.margin = margin(t = 5, b = 5),
          legend.spacing.x = unit(1, "cm"),
          legend.key.size = unit(0.8, "lines"),
          legend.text = element_text(size = 7)) 
dev.off()

```



# Survival analysis with TCGA biolinks

<!-- ```{r} -->

<!-- # Get the events that are significant in the survival analysis for VIPR1 -->
<!-- # survival_assoc_tcga_ids <- -->
<!-- #   te_surv_events %>%  -->
<!-- #   filter(Gene == "AP2A2") %>% -->
<!-- #   pull(Tumor_ids) %>% -->
<!-- #   str_split(",") %>% -->
<!-- #   unlist() %>% -->
<!-- #   unique()  -->


<!-- # Since we don't have AF/AL, do SE events -->

<!-- se_surv_events <- -->
<!-- te_surv_events %>% -->
<!--   filter(event_type == "SE") %>% -->
<!--     # Gene == "VIPR1") %>% -->
<!--   select(Junction, Gene, deltaPSI_TvsN.x, Counts.x, meanPSI_T.x, meanPSI_N.x, deltaPSI_TvsN.x, absDeltaPSI_TvsN.x, starts_with("Survival")) -->

<!-- # Load the complete SE_Cluster data -->
<!-- load(paste0(data_path, "SE_clust_BRCA_TCGA.RData")) -->

<!-- str(clust) -->

<!-- cluster_samples <-  -->
<!--   clust$cluster_exp %>% -->
<!--   filter(Junction %in% se_surv_events$Junction) %>% -->
<!--   filter(Gene == "ECHDC1") %>% -->
<!--   # Get the first three parts of the sample ID -->
<!--   mutate(subject_id = str_split(Sample, "-") %>% map_chr(~paste0(.x[1], "-", .x[2], "-", .x[3]))) %>% -->
<!--   select(subject_id, Classification) -->



<!-- meanPSI_df <- clust$meanPSI -->
<!-- composition_df <- clust$composition -->

<!-- # Convert Cluster values to a common format if necessary -->
<!-- # This step is optional and depends on the actual format of your 'Cluster' column -->

<!-- meanPSI_df$Cluster <- as.factor(meanPSI_df$Cluster) -->
<!-- composition_df$Cluster <- as.numeric(gsub("C", "", composition_df$Cluster)) -->
<!-- composition_df$Cluster <- as.factor(composition_df$Cluster) -->


<!-- # Perform the merge -->
<!-- full_data <- merge(meanPSI_df, composition_df, by = c("Junction", "Gene", "Cluster")) -->



<!-- # Simulate an expanded dataset based on meanPSI and sdPSI -->
<!-- set.seed(123) # For reproducibility -->
<!-- expanded <- do.call(rbind, lapply(1:nrow(clust$meanPSI), function(i) { -->
<!--   data.frame( -->
<!--     Gene = clust$meanPSI$Gene[i], -->
<!--     Cluster = clust$meanPSI$Cluster[i], -->
<!--     meanPSI = rnorm(1000, mean = clust$meanPSI$meanPSI[i], sd = clust$meanPSI$sdPSI[i]) # Generate 1000 points per row as an example -->
<!--   ) -->
<!-- })) -->

<!-- # Plot the expanded dataset -->
<!-- ggplot(expanded, aes(x = meanPSI, y = Gene, fill = as.factor(Cluster))) + -->
<!--   geom_density_ridges() + -->
<!--   theme_bw() + -->
<!--   theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5)) + -->
<!--   labs(x = "Mean PSI", y = "Gene") + # Adjusted label from "Junction" to "Gene" -->
<!--   facet_wrap(~Gene, ncol = 1, scales = "free_y") -->
<!-- dev.off() -->



<!-- # Now do survival analysis with TCGA biolinks -->
<!-- library(TCGAbiolinks) -->

<!-- brca_clinic_gdc <- GDCquery_clinic("TCGA-BRCA", "clinical") -->

<!-- brca_clinic_gdc <- -->
<!-- brca_clinic_gdc %>%  -->
<!--   left_join(cluster_samples, by = c("submitter_id" = "subject_id")) %>% -->
<!--   filter(!is.na(Classification)) -->
<!--   # mutate(has_splice_event = if_else(submitter_id %in% surv_id_splice, "has_event", "lacks_event")) -->

<!-- brca_clinic_gdc %>% -->
<!--   group_by(Classification) %>% -->
<!--   summarise(n = n()) -->

<!-- TCGAanalyze_survival( -->
<!--     data = brca_clinic_gdc, -->
<!--     clusterCol = "Classification", -->
<!--     main = "TCGA Set\n BRCA - ECHDC1", -->
<!--     xlab = "Time since diagnosis (years)", -->
<!--     # risk.table = F, -->
<!--     height = 10, -->
<!--     width=10, -->
<!--     xscale = "d_y" -->
<!-- ) -->






<!-- # In the case of VIPR1, I want to see the combined Junctions -->

<!-- # First, get the Junctions for VIPR1 that have the splicing data  -->
<!-- vipr1_stats <- -->
<!-- te_surv_events %>% -->
<!--   filter(Gene == "VIPR1") %>% -->
<!--   select(Junction, Gene, deltaPSI_TvsN.x, Counts.x, meanPSI_T.x, meanPSI_N.x, deltaPSI_TvsN.x, absDeltaPSI_TvsN.x, starts_with("Survival")) -->




<!-- # DF with Junction, Gene, PSI, Source, Sample, Classification -->
<!-- # Classification is the cluster number from GMM -->
<!-- gmm_class <- AF_clust$cluster_exp %>% as_tibble() -->
<!-- # DF with Junction, Gene, Source, Counts, Percent, Cluster, meanPSI_T, meanPSI_N, deltaPSI_TvsN, absDeltaPSI_TvsN, Wilcox_pval -->
<!-- ts_class <- AF_clust$ts_clusters %>% as_tibble() -->
<!-- composition <- AF_clust$composition %>% as_tibble() -->


<!-- # Check if composition contains VIPR1: -->
<!-- composition %>% -->
<!--   filter(Junction %in% vipr1_stats$Junction) -->

<!-- ts_class %>% -->
<!--   group_by(Junction) %>% -->

<!--   filter(Counts == 463, -->
<!--          Cluster == 1) -->

<!-- # Do a pattern match in composition for Junction that contains chr3:42489477:42489678 -->
<!-- composition %>% -->


<!-- str(AF_clust) -->


<!-- # Get splice junction for AP2A2 -->
<!-- junc <- -->
<!-- te_surv_events %>% -->
<!--   filter(Gene == "AP2A2") %>% -->
<!--   select(Junction) -->


<!-- table(te_surv_events$Junction %in% all_ts_clusters_surv$Junction) -->

<!-- load(paste0(data_path, "AF_clust_BRCA_TCGA.RData")) -->
<!-- # File loaded is AF_clust -->
<!-- AF_clust -->
<!-- AF_250bp_clust -->

<!-- # DF with Junction, Gene, PSI, Source, Sample, Classification -->
<!-- # Classification is the cluster number from GMM -->
<!-- gmm_class <- AF_clust$cluster_exp %>% as_tibble() -->

<!-- # DF with Junction, Gene, Source, Counts, Percent, Cluster, meanPSI_T, meanPSI_N, deltaPSI_TvsN, absDeltaPSI_TvsN, Wilcox_pval -->
<!-- ts_class <- AF_clust$ts_clusters %>% as_tibble() -->

<!-- composition <- AF_clust$composition %>% as_tibble() -->

<!-- composition %>% -->
<!--   filter(Gene == "AP2A2") %>% -->
<!--   distinct(Junction, Gene, Source, Counts, Percent, Cluster) -->
<!--   filter(Junction == junc$Junction) -->

<!-- # Get the Samples fo the classification and the junction -->
<!-- junc_class <-  -->
<!--   ts_class %>% -->
<!--   filter(Junction == junc$Junction) %>% -->
<!--   select(Junction, Classification) %>% -->
<!--   left_join(ts_class, by = "Junction") %>% -->
<!--   select(Junction, Gene, Classification, meanPSI_T, meanPSI_N, deltaPSI_TvsN, absDeltaPSI_TvsN, Wilcox_pval) -->


<!-- gmm_class %>% -->
<!--   filter(Gene == junc$Gene, -->
<!--          PSI > .7 & PSI < .73) %>% -->
<!--   distinct(PSI) -->


<!-- load(paste0(data_path, "all_events_ts_clusters_surv_subpopulation.RData")) -->
<!-- all_ts_clusters_surv %>% -->
<!--   filter(Junction %in% te_surv_events$Junction)  -->

<!-- load(paste0(data_path, "all_events_ts_clusters_surv_250bp_subpopulation.RData")) -->
<!-- all_ts_clusters_surv_250 %>% -->
<!--   filter(Junction %in% te_surv_events$Junction) %>% -->
<!--   filter(Survival_global_pval < 0.05) -->

<!-- signif(as.numeric(junc$absDeltaPSI_TvsN.x), 3) -->

<!--          # Only keep 3 sig figs for deltaPSI -->
<!--          # PSI == signif(as.numeric(junc$absDeltaPSI_TvsN.x), 3)) %>% -->
<!--   distinct(Junction) -->

<!-- gmm_class %>% -->
<!--   filter(PSI == min(PSI)) -->



<!-- library(TCGAbiolinks) -->

<!-- brca_clinic_gdc <- GDCquery_clinic("TCGA-BRCA", "clinical") -->

<!-- # dim(brca_clinic_gdc) -->

<!-- # Submitter id is just first three sections of the barcode -->
<!-- # TCGA-E9-A1NF -->

<!-- surv_id_splice <- -->
<!-- survival_assoc_tcga_ids %>% -->
<!--   str_split("-") %>% -->
<!--   # Select first three sections of the barcode and paste them together -->
<!--   map_chr(~paste0(.x[1], "-", .x[2], "-", .x[3])) -->

<!-- brca_clinic_gdc <- -->
<!-- brca_clinic_gdc %>%  -->
<!--   mutate(has_splice_event = if_else(submitter_id %in% surv_id_splice, "has_event", "lacks_event")) -->

<!-- brca_clinic_gdc %>% -->
<!--   group_by(has_splice_event) %>% -->
<!--   summarise(n = n()) -->

<!-- TCGAanalyze_survival( -->
<!--     data = brca_clinic_gdc, -->
<!--     clusterCol = "has_splice_event", -->
<!--     main = "TCGA Set\n BRCA - VIPR1", -->
<!--     xlab = "Time since diagnosis (years)", -->
<!--     # risk.table = F, -->
<!--     height = 10, -->
<!--     width=10, -->
<!--     xscale = "d_y" -->
<!-- ) -->

<!-- # Do this in a for loop for each gene in te_surv_events -->

<!-- # Get the events that are significant in the survival analysis for VIPR1 -->
<!-- survival_assoc_tcga_ids <- -->
<!--   te_surv_events %>% -->
<!--   pull(Tumor_ids) %>% -->
<!--   str_split(",") %>% -->
<!--   unlist() %>% -->
<!--   unique() -->

<!-- survival_assoc_tcga_ids -->

<!-- # Get the clinical data for BRCA -->
<!-- brca_clinic_gdc <- GDCquery_clinic("TCGA-BRCA", "clinical") -->

<!-- # Submitter id is just first three sections of the barcode -->
<!-- # TCGA-E9-A1NF -->

<!-- surv_id_splice <- -->
<!--   survival_assoc_tcga_ids %>% -->
<!--   str_split("-") %>% -->
<!--   # Select first three sections of the barcode and paste them together -->
<!--   map_chr(~paste0(.x[1], "-", .x[2], "-", .x[3])) -->

<!-- brca_clinic_gdc <- -->
<!--   brca_clinic_gdc %>% -->
<!--   mutate(has_splice_event = if_else(submitter_id %in% surv_id_splice, "has_event", "lacks_event")) -->

<!-- brca_clinic_gdc %>% -->
<!--   group_by(has_splice_event) %>% -->
<!--   summarise(n = n()) -->

<!-- TCGAanalyze_survival( -->
<!--   data = brca_clinic_gdc, -->
<!--   clusterCol = "has_splice_event", -->
<!--   main = "TCGA Set\n BRCA - VIPR1", -->
<!--   height = 10, -->
<!--   width = 10, -->
<!--   filename = "TCGA-BRCA-VIPR1-survival.pdf", -->
<!-- ) -->

<!-- ``` -->


## Make the plots

```{r}



#SKIPPED AND INCLUDED CANCER/NORMAL

pdf(paste0(plots_path, "TE-mediated-splice-class-diverged.pdf"), height = 6, width = 6)

#ready_to_plot <- plottable %>%
#distinct(.keep_all = TRUE) %>% 

ready_to_plot <- all_events_combined %>%
  group_by(event_type) %>% 
    mutate(pos_neg = case_when(TE_IN_REF == "TRUE" & deltaPSI_TvsN.x > 0 ~ "higher in tumor",
                               # TE_IN_ALT == "TRUE" & deltaPSI_TvsN.x < 0 ~ "higher in tumor",
                               TE_IN_REF == "TRUE" & deltaPSI_TvsN.x < 0 ~ "higher in normal",
                               # TE_IN_ALT == "TRUE" & deltaPSI_TvsN.x > 0 ~ "higher in normal"
                               )) %>%
  count(event_type, pos_neg) %>% 
    mutate(plot_val = case_when(pos_neg == "higher in tumor" ~ n,
                                pos_neg == "higher in normal" ~ -n))
  
ready_to_plot %>%
ggplot(aes(x = event_type, y = plot_val, fill = pos_neg)) +
  geom_bar(stat="identity", position="identity")+ 
  geom_text(data = ready_to_plot[ready_to_plot$pos_neg=='higher in normal',], aes(hjust=1.2, label = abs(plot_val))) +
  geom_text(data = ready_to_plot[ready_to_plot$pos_neg=='higher in tumor',], aes(hjust=-.2, label = plot_val)) +
  ggtitle("All quantified TE-mediated splice events") +
  theme_classic() +
  scale_y_continuous(limit = c(-150, 150), 
                     breaks = c(-150, -100, -50, 0, 50, 100, 150),
                     labels = abs(c(-150, -100, -50, 0, 50, 100, 150)))+
  xlab("Splice Event") +
  ylab("Number of events abs(delta PSI > 10%)") +
  theme(axis.text.x = element_text(color = "grey20", size = 12, angle = 0, hjust = .5, vjust = .5, face = "plain"),
        axis.text.y = element_text(color = "grey20", size = 14, angle = 0, hjust = 1, vjust = 0, face = "plain"),  
        axis.title.x = element_text(color = "grey20", size = 14, angle = 0, hjust = .5, vjust = 0, face = "plain"),
        axis.title.y = element_text(color = "grey20", size = 14, angle = 90, hjust = .5, vjust = .5, face = "plain"),
        plot.title = element_text(color = "grey20", size = 16, angle = 0, hjust = .5, vjust = .5, face = "plain")) +
  coord_flip()+
  scale_fill_manual(values=c("#00ADEE","#EC1C24"))
dev.off()

#SKIPPED AND INCLUDED CANCER/NORMAL
pdf(paste0(plots_path, "TE-mediated-splice-repClass-diverged.pdf"), height = 6, width = 6)


ready_to_plot <- all_events_combined %>%
mutate(plot_repClass = (case_when(!(repClass %in% c("DNA", "LTR", "SINE", "LINE")) ~ "Multiple",
                                                                             TRUE ~ repClass))) %>%
  group_by(plot_repClass) %>% 
  mutate(pos_neg = case_when(TE_IN_REF == "TRUE" & deltaPSI_TvsN.x > 0 ~ "higher in tumor",
                             # TE_IN_ALT == "TRUE" & deltaPSI_TvsN.x < 0 ~ "higher in tumor",
                             TE_IN_REF == "TRUE" & deltaPSI_TvsN.x < 0 ~ "higher in normal",
                             # TE_IN_ALT == "TRUE" & deltaPSI_TvsN.x > 0 ~ "higher in normal"
                             )) %>%
  count(plot_repClass, pos_neg) %>% 
  mutate(plot_val = case_when(pos_neg == "higher in tumor" ~ n,
                              pos_neg == "higher in normal" ~ -n))



ready_to_plot %>%
  ggplot(aes(x = plot_repClass, y = plot_val, fill = pos_neg)) +
  geom_bar(stat="identity", position="identity")+ 
  geom_text(data = ready_to_plot[ready_to_plot$pos_neg=='higher in normal',], aes(hjust=1.2, label = abs(plot_val))) +
  geom_text(data = ready_to_plot[ready_to_plot$pos_neg=='higher in tumor',], aes(hjust=-.2, label = plot_val)) +
  scale_y_continuous(limit = c(-200, 200), 
                     breaks = c(-200, -150, -100, -50, 0, 50, 100, 150, 200),
                     labels = abs(c(-200, -150, -100, -50, 0, 50, 100, 150, 200)))+
  ggtitle("All quantified TE-mediated splice events") +
  theme_classic() +
  #scale_fill_brewer(palette="Set1")+
  xlab("Splice Event") +
  ylab("Number of events abs(delta PSI > 10%)") +
  theme(axis.text.x = element_text(color = "grey20", size = 12, angle = 0, hjust = .5, vjust = .5, face = "plain"),
        axis.text.y = element_text(color = "grey20", size = 14, angle = 0, hjust = 1, vjust = 0, face = "plain"),  
        axis.title.x = element_text(color = "grey20", size = 14, angle = 0, hjust = .5, vjust = 0, face = "plain"),
        axis.title.y = element_text(color = "grey20", size = 14, angle = 90, hjust = .5, vjust = .5, face = "plain"),
        plot.title = element_text(color = "grey20", size = 16, angle = 0, hjust = .5, vjust = .5, face = "plain")) +
  coord_flip() +
  scale_fill_manual(values=c("#00ADEE", "#EC1C24"))
dev.off()

#HIGHER IN CANCER

pdf(paste0(plots_path, "TE_higher_cancer_PSI.pdf"), height = 6, width = 6)

total_te_cancer %>%
  mutate(plot_repClass = (case_when(!(repClass %in% c("DNA", "LTR", "SINE", "LINE")) ~ "Multiple",
                                    TRUE ~ repClass))) %>%
  ggplot(aes(x = event_type, fill = plot_repClass)) +
  geom_bar() + ggtitle("TE-mediated splicing events (higher in cancer)") +
  theme_classic() +
  scale_fill_brewer(palette="Set1")+
  xlab("Splice Event") +
  ylab("Number of events (delta PSI > 10%)") +
  theme(axis.text.x = element_text(color = "grey20", size = 12, angle = 0, hjust = .5, vjust = .5, face = "plain"),
        axis.text.y = element_text(color = "grey20", size = 14, angle = 0, hjust = 1, vjust = 0, face = "plain"),  
        axis.title.x = element_text(color = "grey20", size = 14, angle = 0, hjust = .5, vjust = 0, face = "plain"),
        axis.title.y = element_text(color = "grey20", size = 14, angle = 90, hjust = .5, vjust = .5, face = "plain"),
        plot.title = element_text(color = "grey20", size = 16, angle = 0, hjust = .5, vjust = .5, face = "plain")) +
  labs(fill='TE Class')
dev.off()

#HIGHER IN Normal
pdf(paste0(plots_path, "TE_higher_normal_PSI.pdf"), height = 6, width = 6)

total_te_normal %>%
  mutate(plot_repClass = (case_when(!(repClass %in% c("DNA", "LTR", "SINE", "LINE")) ~ "Multiple",
                                    TRUE ~ repClass))) %>%
  ggplot(aes(x = event_type, fill = plot_repClass)) +
  geom_bar() + 
  xlab("Splice Event") +
  ylab("Number of events (delta PSI > 10%)") +
  ggtitle("TE-mediated splice events (higher in normal)")+
  scale_fill_brewer(palette="Set1")+
  theme_classic() +
  #scale_fill_brewer(palette="Set1")+
  theme(axis.text.x = element_text(color = "grey20", size = 12, angle = 0, hjust = .5, vjust = .5, face = "plain"),
        axis.text.y = element_text(color = "grey20", size = 14, angle = 0, hjust = 1, vjust = 0, face = "plain"),  
        axis.title.x = element_text(color = "grey20", size = 14, angle = 0, hjust = .5, vjust = 0, face = "plain"),
        axis.title.y = element_text(color = "grey20", size = 14, angle = 90, hjust = .5, vjust = .5, face = "plain"),
        plot.title = element_text(color = "grey20", size = 18, angle = 0, hjust = .5, vjust = .5, face = "plain")) +
  coord_flip()
dev.off()

#### FIND % of splicing events that overlap TEs

hi_can_status <- total_te_cancer %>% dplyr::select(junction_bak) %>% mutate(TE_status = "Higher_Cancer")
hi_norml_status <- total_te_normal %>% dplyr::select(junction_bak) %>% mutate(TE_status = "Higher_Normal")
hi_status <- rbind(hi_can_status, hi_norml_status)

pdf(paste0(plots_path, "Compare_TE_Status_SUPPA.pdf"), height = 6, width = 6)

left_join(grange_prep, hi_status) %>% 
  mutate(TE_status = case_when(!(TE_status %in% c("Higher_Cancer", "Higher_Normal")) & (deltaPSI_TvsN.x < 0) ~ "No_TE_Higher_Normal",
                               !(TE_status %in% c("Higher_Cancer", "Higher_Normal")) & (deltaPSI_TvsN.x > 0) ~ "No_TE_Higher_Cancer",
                               TRUE ~ TE_status)) %>%
  ggplot(aes(x = event_type, fill = TE_status)) +
  geom_bar() + 
  ggtitle("All SUPPA Quantified events")+
  theme_classic() +
  xlab("Splice Event") +
  ylab("Number of events (delta PSI > 10%)") +
  scale_fill_brewer(palette="Dark2")+
  theme(axis.text.x = element_text(color = "grey20", size = 12, angle = 0, hjust = .5, vjust = .5, face = "plain"),
        axis.text.y = element_text(color = "grey20", size = 14, angle = 0, hjust = 1, vjust = 0, face = "plain"),  
        axis.title.x = element_text(color = "grey20", size = 14, angle = 0, hjust = .5, vjust = 0, face = "plain"),
        axis.title.y = element_text(color = "grey20", size = 14, angle = 90, hjust = .5, vjust = .5, face = "plain"),
        plot.title = element_text(color = "grey20", size = 16, angle = 0, hjust = .5, vjust = .5, face = "plain")) +
  labs(fill='TE Status') +
  coord_flip()
dev.off()


#### TALLY AND PERCENTAGE OF TE-assoc splice events


transcriptome_origin <- left_join(grange_prep, hi_status) %>% 
  mutate(TE_status = case_when(!(TE_status %in% c("Higher_Cancer", "Higher_Normal")) & (deltaPSI_TvsN.x < 0) ~ "No_TE_Higher_Normal",
                               !(TE_status %in% c("Higher_Cancer", "Higher_Normal")) & (deltaPSI_TvsN.x > 0) ~ "No_TE_Higher_Cancer",
                               TRUE ~ TE_status),
         TE_status = factor(TE_status, levels = c("Higher_Cancer", "Higher_Normal", "No_TE_Higher_Normal", "No_TE_Higher_Cancer"))) %>%
  mutate(TE_Overlapping = case_when(TE_status %in% c("Higher_Cancer", "Higher_Normal") ~ "TE overlapping",
                               TRUE ~ "No TE overlapping")) %>%
  distinct(junction_bak, .keep_all = T) %>%
  group_by(TE_status, TE_Overlapping, transcriptome_origin) %>%
  tally(name = "count") %>%
  mutate(transcriptome_origin, percent = (count/sum(count))*100) %>%
  mutate(cancer_normal_expression = case_when(TE_status == "Higher_Cancer" ~ "Higher Cancer",
                               TE_status == "Higher_Normal" ~ "Higher Normal",
                               TE_status == "No_TE_Higher_Normal" ~ "Higher Normal",
                               TE_status == "No_TE_Higher_Cancer" ~ "Higher Cancer",
                               TRUE ~ "No_TE")) %>%
  mutate(transcriptome_origin = factor(transcriptome_origin, c("LR-seq", "GENCODE")))

pdf(paste0(plots_path, "Transcriptome_origin.pdf"), height = 6, width = 6)

transcriptome_origin %>%
  ggplot(aes(x = cancer_normal_expression, y = count, fill = transcriptome_origin, alpha = TE_Overlapping)) +
  geom_bar(stat = "identity", position = position_dodge(width = 0.9)) + 
  geom_text(aes(label = count, y = count + 70),
            position = position_dodge(width = 0.9), size = 4) +

  ggtitle("All differential splicing events")+
  coord_flip()+
  theme_classic() +
  xlab("") +
  ylab("Number of events abs(delta PSI > 10%)") +
  scale_fill_manual(values = c("LR-seq" = "#68457A", "GENCODE" = "#FF9933")) +
  scale_alpha_manual(values = c("No TE overlapping" = 0.5, "TE overlapping" = 1)) +
  theme(axis.text.x = element_text(color = "grey20", size = 12, angle = 0, hjust = .5, vjust = .5, face = "plain"),
        axis.text.y = element_text(color = "grey20", size = 14, angle = 0, hjust = 1, vjust = 0, face = "plain"),  
        axis.title.x = element_text(color = "grey20", size = 14, angle = 0, hjust = .5, vjust = 0, face = "plain"),
        axis.title.y = element_text(color = "grey20", size = 14, angle = 90, hjust = .5, vjust = .5, face = "plain"),
        plot.title = element_text(color = "grey20", size = 16, angle = 0, hjust = .5, vjust = .5, face = "plain")) +
  labs(fill='Transcriptome') +
  coord_flip() +
  facet_wrap(~TE_Overlapping, ncol = 1, scales = "free_y")
dev.off()



# Do a pathway enrichment for the TE higher in cancer or te higher in normal events


library(enrichR)

alldbs <- listEnrichrDbs()

dbs <- c("KEGG_2019_Human", "Reactome_2016", "GO_Biological_Process_2018", "GO_Cellular_Component_2018", "GO_Molecular_Function_2018", "MSigDB_Hallmark_2020", "MSigDB_Oncogenic_Signatures")

# Get the gene names for the TE higher in cancer and TE higher in normal events
te_cancer_genes <- total_te_cancer %>% pull(Gene) %>% unique()
te_normal_genes <- total_te_normal %>% pull(Gene) %>% unique()

enriched_cancer <- enrichr(te_cancer_genes, dbs)
enriched_normal <- enrichr(te_normal_genes, dbs)

plotEnrich(enriched_cancer[[7]], showTerms = 20, numChar = 40, y = "Count", orderBy = "P.value")
plotEnrich(enriched_normal[[7]], showTerms = 20, numChar = 40, y = "Count", orderBy = "P.value")


plotEnrich(enriched_cancer[[6]], showTerms = 20, numChar = 40, y = "Count", orderBy = "P.value")
plotEnrich(enriched_normal[[6]], showTerms = 20, numChar = 40, y = "Count", orderBy = "P.value")


# Make a venn diagram of the genes
library(ggvenn)
venn_genes <- list(cancer_genes = te_cancer_genes, normal_genes = te_normal_genes)
venn_genes_plot <- ggvenn(venn_genes)


library(patchwork)

venn_genes_plot





pdf(paste0(plots_path, "Transcriptome_origin_wide.pdf"), height = 3, width = 6)

transcriptome_origin %>%
  mutate(cancer_normal_expression = case_when(TE_status == "Higher_Cancer" ~ "Cancer",
                               TE_status == "Higher_Normal" ~ "Normal",
                               TE_status == "No_TE_Higher_Normal" ~ "Normal",
                               TE_status == "No_TE_Higher_Cancer" ~ "Cancer",
                               TRUE ~ "No_TE")) %>%
  ggplot(aes(x = cancer_normal_expression, y = count, fill = transcriptome_origin, alpha = TE_Overlapping)) +
  geom_bar(stat = "identity", position = position_dodge(width = 0.9)) + 
  geom_text(aes(label = count, y = count + 70),
            position = position_dodge(width = 0.9), size = 4) +
  ggtitle("All differential splicing events")+
  # coord_flip()+
  theme_classic() +
  xlab("") +
  ylab("Number of events abs(delta PSI > 10%)") +
  scale_fill_manual(values = c("LR-seq" = "#68457A", "GENCODE" = "#FF9933")) +
  scale_alpha_manual(values = c("No TE overlapping" = 0.5, "TE overlapping" = 1)) +
  theme(axis.text.x = element_text(color = "grey20", size = 12, angle = 0, hjust = .5, vjust = .5, face = "plain"),
        axis.text.y = element_text(color = "grey20", size = 14, angle = 0, hjust = 1, vjust = 0, face = "plain"),  
        axis.title.x = element_text(color = "grey20", size = 14, angle = 0, hjust = .5, vjust = 0, face = "plain"),
        axis.title.y = element_text(color = "grey20", size = 14, angle = 90, hjust = .5, vjust = .5, face = "plain"),
        plot.title = element_text(color = "grey20", size = 16, angle = 0, hjust = .5, vjust = .5, face = "plain")) +
  labs(fill='Transcriptome') +
  # coord_flip() +
  facet_wrap(~TE_Overlapping, ncol = 2, scales = "free_y")

dev.off()


```


## Check genes with TE-associated splice events
```{r}
annotated_set <-
left_join(grange_prep, hi_status) %>% 
  mutate(TE_status = case_when(!(TE_status %in% c("Higher_Cancer", "Higher_Normal")) & (deltaPSI_TvsN.x < 0) ~ "No_TE_Higher_Normal",
                               !(TE_status %in% c("Higher_Cancer", "Higher_Normal")) & (deltaPSI_TvsN.x > 0) ~ "No_TE_Higher_Cancer",
                               TRUE ~ TE_status),
         TE_status = factor(TE_status, levels = c("Higher_Cancer", "Higher_Normal", "No_TE_Higher_Normal", "No_TE_Higher_Cancer"))) %>%
  mutate(TE_Overlapping = case_when(TE_status %in% c("Higher_Cancer", "Higher_Normal") ~ "TE overlapping",
                               TRUE ~ "No TE overlapping")) %>%
  distinct(junction_bak, .keep_all = T) %>%
  filter(TE_Overlapping == "TE overlapping")


# get the unique genes for each (Higher Cancer or Higher Normal) - Then check the overlap with the other set


# Make a 4 way venn diagram

annotated_categories <-
left_join(grange_prep, hi_status) %>% 
  mutate(TE_status = case_when(!(TE_status %in% c("Higher_Cancer", "Higher_Normal")) & (deltaPSI_TvsN.x < 0) ~ "No_TE_Higher_Normal",
                               !(TE_status %in% c("Higher_Cancer", "Higher_Normal")) & (deltaPSI_TvsN.x > 0) ~ "No_TE_Higher_Cancer",
                               TRUE ~ TE_status),
         TE_status = factor(TE_status, levels = c("Higher_Cancer", "Higher_Normal", "No_TE_Higher_Normal", "No_TE_Higher_Cancer"))) %>%
  mutate(TE_Overlapping = case_when(TE_status %in% c("Higher_Cancer", "Higher_Normal") ~ "TE overlapping",
                               TRUE ~ "No TE overlapping")) %>%
  distinct(junction_bak, .keep_all = T)

unique_genes <-
annotated_categories %>%
  dplyr::select(Gene, TE_status) %>%
  group_by(TE_status) %>%
  distinct(Gene)

# Convert the dataframe into a list of 4 categories based on TE_status
unique_genes <- unique_genes %>%
  split(.$TE_status) %>%
  lapply(function(x) x$Gene)

ggVennDiagram(unique_genes)

# make an upset plot
library(UpSetR)

pdf(paste0(plots_path, "Upset_plot.pdf"), height = 6, width = 6)

# Your upset function call with color settings
# Assuming the color codes for LR-seq and GENCODE are as follows
color1 <- "#6C3483" # A purple color similar to LR-seq
color2 <- "#F39C12" # An orange color similar to GENCODE

library(ComplexHeatmap)

m <- make_comb_mat(unique_genes)

m <- m[comb_degree(m) %in% c(1,2)]

pdf(paste0(plots_path, "Upset_plot.pdf"), height = 3, width = 6)

UpSet(m, 
      set_order = order(set_size(m), decreasing = TRUE),
      comb_order = order(comb_size(m), decreasing = T),
      
      )
dev.off()

?UpSet

    # comb_col = c("red", "blue", "black")[comb_degree(m)])
# 
# ?UpSet
# 
# UpSet(unique_genes, sets = c("Higher_Cancer", "Higher_Normal", "No_TE_Higher_Normal", "No_TE_Higher_Cancer"), 
#       sets.bar.color = c(color1, color1, color2, color2), order.by = "freq", main.bar.color = "black", matrix.color = "black")


unique_genes <- annotated_set %>%
  distinct(TE_status, Gene)

library(ggVennDiagram)

pdf(paste0(plots_path, "Unique_genes.pdf"), height = 6, width = 6)

# Make a list of vectors using each group of TE_status
unique_genes %>%
  mutate(TE_status = as.character(TE_status)) %>%
    # mutate(TE_status = case_when(TE_status == "Higher_Normal" ~ "TE PSI higher in normal"),
    #      TE_status = case_when(TE_status == "Higher_Cancer" ~ "TE PSI higher in cancer")) %>%
  group_by(TE_status) %>%
  split(.$TE_status) %>%
  lapply(function(x) x$Gene) %>%
  ggVennDiagram(label = "both" ) +
  coord_flip()
# ?ggVennDiagram
shared_genes <- annotated_set %>%
  distinct(Gene, TE_status) %>%
  filter(TE_status %in% c("Higher_Cancer", "Higher_Normal")) %>%
  group_by(Gene) %>%
  mutate(count = dplyr::n()) %>%
  dplyr::filter(count > 1) %>%
  distinct(Gene, .keep_all = T)

#now get the unique genes
unique_genes <- annotated_set %>%
  distinct(Gene, TE_status) %>%
  filter(!(Gene %in% shared_genes$Gene))

# Now do an enrichment using enrichR in each category
# https://cran.r-project.org/web/packages/enrichR/vignettes/enrichR.html

library(enrichR)

# Get the gene symbols
gene_symbols_higher_norma <- unique_genes %>%
  filter(TE_status == "Higher_Normal") %>%
  distinct(Gene) %>%
  pull(Gene)

gene_symbols_higher_cancer <- unique_genes %>%
  filter(TE_status == "Higher_Cancer") %>%
  distinct(Gene) %>%
  pull(Gene)

gene_symbols_shared <-
  shared_genes %>%
  distinct(Gene) %>%
  pull(Gene)

# Now do the enrichment

dbs <- c("GO_Molecular_Function_2015", "GO_Cellular_Component_2015", "GO_Biological_Process_2021", "MSigDB_Oncogenic_Signatures")

enriched_normal <- enrichr(gene_symbols_higher_norma, dbs)
enriched_cancer <- enrichr(gene_symbols_higher_cancer, dbs)
enriched_shared <- enrichr(gene_symbols_shared, dbs)

# Plot the results for each category
plotEnrich(enriched_normal[[4]], 
           showTerms = 20, 
           numChar = 40, 
           y = "Count", 
           orderBy = "P.value")

plotEnrich(enriched_cancer[[4]],
           showTerms = 20, 
           numChar = 40, 
           y = "Count", 
           orderBy = "P.value")

plotEnrich(enriched_shared[[4]],
           showTerms = 20, 
           numChar = 40, 
           y = "Count", 
           orderBy = "P.value")


```


##Splicing quantifications that overlap TEs

## Load packages and data
```{r}

library(tidyverse)
library(ggplot2)
library(cowplot)
library(ggplotify)
library(pheatmap)
library(RColorBrewer)
library(rtracklayer)
library(plyranges)

## Load TEs higher in cancer
cancer_hits <- read_tsv(paste0(tables_path, "TE_higher_cancer.tsv")) %>%
  mutate("repOrientation" = (case_when(
                                    strand == repStrand ~ "sense",
                                    strand != repStrand ~ "antisense")))

gencode_first <- import.gff2(paste0(data_path, "first_exons_gencode_TSL1.gtf"))
gencode_last <- import.gff2(paste0(data_path, "last_exons_gencode_TSL1.gtf"))

```


#Functions for AF hits
```{r}

total_AF <- cancer_hits %>% 
  distinct(Junction, .keep_all = TRUE) %>%
  filter(number_controls_all > 10,
         Counts.x > 30) %>%
  filter(event_type == "AF")


# FOR THE POSITIVE DPSI EVENTS
AF_grange_second_jxn <- total_AF %>% 
                        filter(deltaPSI_TvsN.x > 0) %>%
                        dplyr::select(-start, -end) %>% 
                        makeGRangesFromDataFrame(seqnames.field = "seqnames",
                                                 start.field = "start_2",
                                                 end.field = "end_2",
                                                keep.extra.columns = T)

## FOR THE NEGATIVE DPSI EVENTS 
AF_grange_first_jxn <- total_AF %>% 
  filter(deltaPSI_TvsN.x < 0) %>%
  makeGRangesFromDataFrame(seqnames.field = "seqnames",
                           keep.extra.columns = T)


twofifty_5p_gencode_overlap <- function(granges_object){
  ### GET 5p and 3p for the events
  granges_object_5p <- mutate(anchor_5p(granges_object), width = 1)
  granges_object_3p <- mutate(anchor_3p(granges_object), width = 1)
  
  ### GET 5p and 3p for GENCODE FIRST exons
  gencode_first_5p <- mutate(anchor_5p(gencode_first), width = 1)
  gencode_first_3p <- mutate(anchor_3p(gencode_first), width = 1)
  
  ### INTERSECT 3p end maxgap == 2
  fiveprime_overlaps <- find_overlaps(granges_object_3p, gencode_first_3p, maxgap = 2)
  
  ### INTERSECT 5p end maxgap == 250
  threeprime_overlaps <- find_overlaps(granges_object_5p, gencode_first_5p, maxgap = 250)

  final_grange <- granges_object %>% filter(junction_bak %in% fiveprime_overlaps$junction_bak & junction_bak %in% threeprime_overlaps$junction_bak)
  
  return(final_grange)
}

twofifty_3p_gencode_overlap <- function(granges_object){
  ### GET 5p and 3p for the events
  granges_object_5p <- mutate(anchor_5p(granges_object), width = 1)
  granges_object_3p <- mutate(anchor_3p(granges_object), width = 1)
  
  ### GET 5p and 3p for GENCODE FIRST exons
  gencode_last_5p <- mutate(anchor_5p(gencode_last), width = 1)
  gencode_last_3p <- mutate(anchor_3p(gencode_last), width = 1)
  
  ### INTERSECT 3p end maxgap == 2
  fiveprime_overlaps <- find_overlaps(granges_object_3p, gencode_last_3p, maxgap = 250)
  
  ### INTERSECT 5p end maxgap == 250
  threeprime_overlaps <- find_overlaps(granges_object_5p, gencode_last_5p, maxgap = 2)
  
  final_grange <- granges_object %>% filter(junction_bak %in% fiveprime_overlaps$junction_bak & junction_bak %in% threeprime_overlaps$junction_bak)
  
  return(final_grange)
}

```

## First set of TE overlaps
```{r}


## GET OVERLAPPING HITS
result_AF_pos <- twofifty_5p_gencode_overlap(AF_grange_second_jxn)
result_AF_neg <- twofifty_5p_gencode_overlap(AF_grange_first_jxn)
AF_pos_neg_res <- bind_ranges(result_AF_pos, result_AF_neg)

## GET OVERLAPPING from original GRANGES
all_TE_af_vs_gencode <-total_AF %>% filter(junction_bak %in% AF_pos_neg_res$junction_bak)



## RI Filter
extract_unique_classes <- function(repClass) {
  # Split the string by comma
  split_classes <- strsplit(repClass, ",")
  
  # Get unique classes for each element
  unique_classes <- lapply(split_classes, function(x) {
    # Remove any leading or trailing whitespace
    x <- trimws(x)
    
    # Get unique classes
    unique_x <- unique(x)
    
    # Sort the unique classes alphabetically
    unique_x <- sort(unique_x)
    
    # Collapse back into a single string
    paste(unique_x, collapse = ",")
  })
  
  # Return the result as a vector
  unlist(unique_classes)
}




total_SE <- cancer_hits %>% 
  distinct(Junction, .keep_all = TRUE) %>%
  filter(number_controls_all > 10,
         Counts.x > 30) %>%
  filter(event_type == "SE")

total_RI <- cancer_hits %>% 
  distinct(Junction, .keep_all = TRUE) %>%
  filter(number_controls_all > 10,
         Counts.x > 30) %>%
  filter(event_type == "RI") %>%
  mutate(repClass = extract_unique_classes(repClass))
  # we need to collapse the classes from the list so that there can be only one of each

total_AL <- cancer_hits %>% 
  distinct(Junction, .keep_all = TRUE) %>%
  filter(number_controls_all > 10,
         Counts.x > 30) %>%
  filter(event_type == "AL")

# FOR THE POSITIVE DPSI EVENTS
AL_grange_second_jxn <- total_AL %>% 
  filter(deltaPSI_TvsN.x > 0) %>%
  dplyr::select(-start, -end) %>% 
  makeGRangesFromDataFrame(seqnames.field = "seqnames",
                           start.field = "start_2",
                           end.field = "end_2",
                           keep.extra.columns = T)

## FOR THE NEGATIVE DPSI EVENTS 
AL_grange_first_jxn <- total_AL %>% 
  filter(deltaPSI_TvsN.x < 0) %>%
  makeGRangesFromDataFrame(seqnames.field = "seqnames",
                           keep.extra.columns = T)


result_AL_pos <- twofifty_3p_gencode_overlap(AL_grange_second_jxn)
result_AL_neg <- twofifty_3p_gencode_overlap(AL_grange_first_jxn)


AL_pos_neg_res <- bind_ranges(result_AL_pos, result_AL_neg)

## GET OVERLAPPING from original GRANGES
all_TE_al_vs_gencode <-total_AL %>% filter(junction_bak %in% AL_pos_neg_res$junction_bak)

```

## AF and AL plots
``` {r}

x <- total_AF %>%
  group_by(repClass, repOrientation) %>%
  tally() %>%
  ggplot(aes(x = repClass, y=n, fill = repOrientation)) +
  geom_bar(stat = "identity") +
  geom_text(
    aes(label = stat(y), group = repClass), 
    stat = 'summary', fun = sum, vjust = -1) +
  theme_classic() +
  scale_fill_manual(values=c("sense" = "#606060", "antisense"="#D8E1E7")) +
  ylab("Number of splice events")+
  xlab("TE class") +
  labs(fill="TE orientation\nvs. gene") +
  ggtitle("AF events")+
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

y <- total_SE %>%
  group_by(repClass, repOrientation) %>%
  tally() %>%
  ggplot(aes(x = repClass, y=n, fill = repOrientation)) +
  geom_bar(stat = "identity") +
  geom_text(
    aes(label = stat(y), group = repClass), 
    stat = 'summary', fun = sum, vjust = -1) +
  theme_classic() +
  scale_fill_manual(values=c("sense" = "#606060", "antisense"="#D8E1E7")) +
  ylab("Number of splice events")+
  xlab("TE class") +
  labs(fill="TE orientation\nvs. gene") +
  ggtitle("CA events")+
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

z <- total_AL %>%
  group_by(repClass, repOrientation) %>%
  tally() %>%
  ggplot(aes(x = repClass, y=n, fill = repOrientation)) +
  geom_bar(stat = "identity") +
  geom_text(
    aes(label = stat(y), group = repClass), 
    stat = 'summary', fun = sum, vjust = -1) +
  theme_classic() +
  scale_fill_manual(values=c("sense" = "#606060", "antisense"="#D8E1E7")) +
  ylab("Number of splice events")+
  xlab("TE class") +
  labs(fill="TE orientation\nvs. gene") +
  ggtitle("AL events")+
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

zz <- 
  total_RI %>%
  mutate(repClass = if_else(str_detect(repClass, ","), "Multiple\nRepeats", repClass),
         repOrientation = if_else(
                                  str_detect(repStrand, ","),
                                  "Multiple\nOrientations",
                                  repOrientation),
         repOrientation = if_else(
           is.na(repOrientation), 
           "Multiple\nOrientations", 
           repOrientation
         )) %>%
  group_by(repClass, repOrientation) %>%
  tally() %>%
  ggplot(aes(x = repClass, y=n, fill = repOrientation)) +
  geom_bar(stat = "identity") +
  geom_text(
    aes(label = stat(y), group = repClass), 
    stat = 'summary', fun = sum, vjust = -1) +
  theme_classic() +
  scale_fill_manual(values=c("sense" = "#606060", 
                             "antisense"="#D8E1E7", 
                             "Multiple\nOrientations" = "black")) +
  ylab("Number of splice events")+
  xlab("TE class") +
  labs(fill="TE orientation\nvs. gene") +
  ggtitle("RI events") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))+
  theme(legend.position = "bbottom")


```


## AF TE switches
```{r}


top_each_AF <- all_TE_af_vs_gencode %>%  
  as_tibble() %>%
  arrange(desc(Counts.x)) %>%
  distinct(gene, .keep_all = TRUE) %>%
  mutate_at(vars(contains("Gene")), (as.factor)) %>% 
  filter(number_controls_all > 10,
         Counts.x > 30) %>%
  slice(1:5) %>%
  mutate(x_order = seq(n()))



AF_FILTERED_EVENTS <- cancer_hits %>% 
  distinct(gene, .keep_all = TRUE) %>%
  arrange(desc(Counts.x)) %>%
  mutate_at(vars(contains("Gene")), (as.factor)) %>% 
  filter(number_controls_all > 10,
         Counts.x > 30,
         absDeltaPSI_TvsN.x >= .1) %>%
  filter(event_type == "AF") 


```


## Plot all of the isoform switches
```{r}

a <- 
top_each_AF %>%
  ungroup() %>%
  mutate(Counts.x = as.numeric(Counts.x)) %>%
  mutate(gene_te = paste0(Gene, " - ", repClass)) %>%
  mutate(gene_te = forcats::fct_reorder(gene_te, Counts.x)) %>%
  ggplot(aes(x = gene_te, y = Counts.x, fill = deltaPSI_TvsN.x)) +
  geom_bar(stat = "identity") +
  scale_fill_gradient2(
    limits = c(-1.0, 1.0),
    breaks = c(1.0, 0.5, 0, -0.5, -1.0),
    low = "#00478F",
    mid = "white",
    high = "#FF5D00",
    midpoint = 0,
    na.value = NA
  ) +
  coord_flip() +
  theme_classic() +
  ylab("Total patients with splice event") +
  xlab("Gene symbol - TE class") +
  labs(fill = "Percent spliced in\n +Tumor/-Normal") +
  ggtitle("AF events")+
  theme(legend.position = "none")



top_each_SE <- cancer_hits %>% 
  distinct(gene, .keep_all = TRUE) %>%
  arrange(desc(Counts.x)) %>%
  mutate_at(vars(contains("Gene")), (as.factor)) %>% 
  filter(number_controls_all > 10,
         Counts.x > 30,
         deltaPSI_TvsN.x >= .1) %>%
  filter(event_type == "SE") %>%
  slice(1:5)%>%
  mutate(x_order = seq(n()))

top_each_SE %>% dplyr::select(seqnames, start, end, Gene, start_2, end_2)


b <- top_each_SE %>% 
  ungroup() %>%
  mutate(Counts.x = as.numeric(Counts.x)) %>%
  mutate(gene_te = paste0(Gene, " - ", repClass)) %>%
  mutate(gene_te = forcats::fct_reorder(gene_te, Counts.x)) %>%
  ggplot(aes(x = gene_te, y = Counts.x, fill = deltaPSI_TvsN.x)) +
  geom_bar(stat = "identity") +
  scale_fill_gradient2(
    limits = c(-1.0, 1.0),
    breaks = c(1.0, 0.5, 0, -0.5, -1.0),
    low = "#00478F",
    mid = "white",
    high = "#FF5D00",
    midpoint = 0,
    na.value = NA
  ) +  coord_flip() +
  theme_classic() +
  ylab("Total patients with splice event")+
  xlab("Gene symbol - TE class") +
  labs(fill="Percent spliced in\n +Tumor/-Normal")+
  ggtitle("CA events")+
  theme(legend.position = "none")

top_each_AL <- all_TE_al_vs_gencode %>%  
  as_tibble() %>%
  arrange(desc(Counts.x)) %>%
  distinct(gene, .keep_all = TRUE) %>%
  mutate_at(vars(contains("Gene")), (as.factor)) %>% 
  filter(number_controls_all > 10,
         Counts.x > 30) %>%
  slice(1:5) %>%
  mutate(x_order = seq(n()))


AL_FILTERED_EVENTS <- cancer_hits %>% 
  distinct(gene, .keep_all = TRUE) %>%
  arrange(desc(Counts.x)) %>%
  mutate_at(vars(contains("Gene")), (as.factor)) %>% 
  filter(number_controls_all > 10,
         Counts.x > 30,
         absDeltaPSI_TvsN.x > .1) %>%
  filter(event_type == "AL")

#AL_FILTERED_EVENTS %>% write_tsv("AL_FILTERED_EVENTS.tsv")
top_each_AL %>% dplyr::select(seqnames, start, end, Gene, start_2, end_2)

c <- top_each_AL %>% 
  ungroup() %>%
  mutate(Counts.x = as.numeric(Counts.x)) %>%
  mutate(gene_te = paste0(Gene, " - ", repClass)) %>%
  mutate(gene_te = forcats::fct_reorder(gene_te, Counts.x)) %>%
  ggplot(aes(x = gene_te, y = Counts.x, fill = deltaPSI_TvsN.x)) +
  geom_bar(stat = "identity") +
  scale_fill_gradient2(
    limits = c(-1.0, 1.0),
    breaks = c(1.0, 0.5, 0, -0.5, -1.0),
    low = "#00478F",
    mid = "white",
    high = "#FF5D00",
    midpoint = 0,
    na.value = NA
  ) +  coord_flip() +
  theme_classic() +
  ylab("Total patients with splice event")+
  xlab("Gene symbol - TE class") +
  labs(fill="Percent spliced in\n +Tumor/-Normal")+
  ggtitle("AL events") +
  theme(legend.position = "none")

RI_FILTERED_EVENTS <- cancer_hits %>% 
  distinct(gene, .keep_all = TRUE) %>%
  arrange(desc(Counts.x)) %>%
  mutate_at(vars(contains("Gene")), (as.factor)) %>% 
  filter(number_controls_all > 10,
         Counts.x > 30,
         absDeltaPSI_TvsN.x > .1) %>%
  filter(event_type == "RI")

d <- 
  total_RI %>% 
  ungroup() %>%
    mutate(repClass = if_else(str_detect(repClass, ","), "Multiple\nRepeats", repClass),
         repOrientation = if_else(
                                  str_detect(repStrand, ","),
                                  "Multiple\nOrientations",
                                  repOrientation),
         repOrientation = if_else(
           is.na(repOrientation), 
           "Multiple\nOrientations", 
           repOrientation
         )) %>%
  mutate(Counts.x = as.numeric(Counts.x)) %>%
  mutate(gene_te = paste0(Gene, " - ", repClass)) %>%
  mutate(gene_te = forcats::fct_reorder(gene_te, Counts.x)) %>%
  arrange(desc(gene_te))%>%
    dplyr::slice(1:5) %>%
  ggplot(aes(x = gene_te, y = Counts.x, fill = deltaPSI_TvsN.x)) +
  geom_bar(stat = "identity") +
  scale_fill_gradient2(
    limits = c(-1.0, 1.0),
    breaks = c(1.0, 0.5, 0, -0.5, -1.0),
    low = "#00478F",
    mid = "white",
    high = "#FF5D00",
    midpoint = 0,
    na.value = NA
  ) +  coord_flip() +
  theme_classic() +
  ylab("Total patients with splice event")+
  xlab("Gene symbol - TE class") +
  labs(fill="Percent spliced in\n +Tumor/-Normal")+
  ggtitle("RI events") +
  theme(legend.position = "bottom")



pdf(paste0(plots_path, "poster_top_ALL.pdf"), height = 5.25, width = 8)

# plot_grid(x,y,z, a,b,c,  ncol = 3, nrow = 2)

# Use patchwork to combine the plots
# Also combine the legends
(x | y | z | zz) /
(a | b | c | d) &
  theme(legend.position = "none",
        axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        axis.text.x = element_text(angle = 45, hjust = 1),
        axis.text.y = element_text(angle = -45, vjust = 1))

# ) /
#   (a + b + c + d &theme(legend.position = "none",
#         axis.title.x = element_blank(),
#         axis.title.y = element_blank())) &
#   plot_layout(byrow = TRUE,
#               guides = "collect",  # Collect the legends
#               axes = "collect",
#               tag_level = "new",
#               axis_titles = "collect")


dev.off()




hits_to_validate <- bind_rows(top_each_AF, top_each_SE, top_each_AL)
hits_to_validate %>% write_tsv(paste0(tables_path, "top_splicing_hits_to_validate.tsv"))

```


## Stats tests
```{r}
## Statement to prove - Alternative cassette exons are the result of the polyT branch site in SINEs and LINEs and their orientation with respect to the gene TSS. This is overrepresented compared to LTR and DNA TEs.

## Null hypothesis - Alternative cassette exons are not the result of the polyT branch site in SINEs and LINEs, and are just the result of the higher number of SINES and LINEs in the introns of genes.

## Step 1: Calculate the number of SINEs, LINES, DNA, and LTRs between a TSS and TTS of a gene in gencode


# Load gencode
import.gff2(paste0(data_path, "gencode.sorted.gtf")) -> gencode

#select protein coding genes in gencode
gencode_genes <-
  gencode %>% 
    filter(type == "gene",
           gene_type == "protein_coding") %>%
  mutate(gene_length = end - start)


#intersect the simp_rmsk with the intronic regions of genes
simp_rmsk <-
simp_rmsk %>%
  mutate(repStart = start(simp_rmsk),
         repEnd = end(simp_rmsk))

find_overlaps(gencode_genes, simp_rmsk, suffix = c(".gencode", ".tes")) -> te_gene_overlaps

## Step 2: Determine the TE strand orientation with respesct to the gene TSS

te_gene_overlaps_oriented <-
te_gene_overlaps %>% 
  as_tibble() %>%
  mutate(te_orientation = case_when(strand == repStrand ~ "sense",
                                    strand != repStrand ~ "antisense",
                                    T ~ "NA")) %>%
  dplyr::select(seqnames, start, end, width, gene_length, strand, source, gene_id, gene_type, gene_name, level, repStart, repEnd, repName, repClass, repFamily, repStrand, te_orientation) 

## Count the number of TEs in each orientation

count_of_tes_oriented <-
te_gene_overlaps_oriented %>%
  group_by(gene_id, gene_length, te_orientation, repClass) %>%
  count(name = "number_of_tes") %>%
  mutate(normalized_number_of_tes = number_of_tes/gene_length)


## Step 3: Get diff splicing events

included_cassette_te_events <-
total_SE %>%
  group_by(repClass, repOrientation) %>%
  tally()


total_SE %>%
  filter(repClass == "SINE",
         repOrientation == "antisense") %>%
  group_by(repFamily) %>%
  tally()

## Step 4: contingency table

# Aggregate counts from count_of_tes_oriented dataset
tes_totals <- count_of_tes_oriented %>%
  group_by(repClass, te_orientation) %>%
  summarize(total_tes = sum(number_of_tes)) %>%
  ungroup()



# Create an empty 4x4 matrix for the contingency table
contingency_table <- matrix(0, nrow=4, ncol=4, 
                            dimnames = list(unique(tes_totals$repClass), 
                                            c("antisense_all_genes", "sense_all_genes", "antisense_splice_included", "sense_splice_included")))

# Fill the table with values from tes_totals
for (i in 1:nrow(tes_totals)) {
  column_name <- paste(tes_totals$te_orientation[i], "all_genes", sep="_")
  contingency_table[tes_totals$repClass[i], column_name] <- tes_totals$total_tes[i]
}

# Add values from included_cassette_te_events to the table
for (i in 1:nrow(included_cassette_te_events)) {
  column_name <- paste(included_cassette_te_events$repOrientation[i], "splice_included", sep="_")
  contingency_table[included_cassette_te_events$repClass[i], column_name] <- included_cassette_te_events$n[i]
}

contingency_table

# Conduct the chi-squared test on the contingency table
chi2_test_result <- chisq.test(contingency_table)
# Convert the contingency table to a dataframe and retain rownames
df_contingency <- as.data.frame(contingency_table) %>%
  rownames_to_column(var = "repClass")

# Convert the contingency table to a long format suitable for ggplot2
long_data <- df_contingency %>%
  pivot_longer(cols = -repClass, names_to = "Type", values_to = "Count")


# Plot
ggplot(long_data, aes(x = Type, y = Count, fill = repClass)) +
  geom_bar(stat = "identity", position = "stack") +
  facet_wrap(~repClass, scales = "free_y") +
  theme_minimal() +
  labs(title = "Distribution of repClass TEs based on Orientation", 
       x = "Orientation Type", 
       y = "Count of TEs", 
       fill = "RepClass") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))






####

ratio_per_gene <- count_of_tes_oriented %>%
  group_by(gene_id, repClass) %>%
  summarize(antisense_count = sum(ifelse(te_orientation == "antisense", number_of_tes, 0)),
            sense_count = sum(ifelse(te_orientation == "sense", number_of_tes, 0))) %>%
  mutate(gene_ratio = antisense_count / sense_count) %>%
  filter(sense_count != 0)  # Avoid dividing by zero

average_expected_ratios <- ratio_per_gene %>%
  group_by(repClass) %>%
  summarize(expected_ratio = mean(gene_ratio, na.rm = TRUE))

print(average_expected_ratios)



observed_ratios <- included_cassette_te_events %>%
  group_by(repClass) %>%
  summarize(antisense_count = sum(ifelse(repOrientation == "antisense", n, 0)),
            sense_count = sum(ifelse(repOrientation == "sense", n, 0))) %>%
  mutate(observed_ratio = antisense_count / sense_count) %>%
  filter(sense_count != 0)  # Avoid dividing by zero

print(observed_ratios)


combined_ratios <- left_join(observed_ratios[,c('repClass', 'observed_ratio')],
                             average_expected_ratios, 
                             by = "repClass")

print(combined_ratios)

#Calculate a p value from the expected and observed ratios
combined_ratios <- combined_ratios %>%
  mutate(p_value = ifelse(observed_ratio > expected_ratio, 
                          pchisq(observed_ratio, df = 1, lower.tail = FALSE), 
                          pchisq(observed_ratio, df = 1, lower.tail = TRUE)))


library(ggplot2)

# Convert data to long format
long_ratios <- combined_ratios %>%
  pivot_longer(cols = c(observed_ratio, expected_ratio), 
               names_to = "Ratio_Type", 
               values_to = "Value")


# Plot
ggplot(long_ratios, aes(x = repClass, y = Value, fill = Ratio_Type)) +
  geom_bar(stat = "identity", position = "dodge") +
  theme_minimal() +
  labs(title = "Comparison of Observed vs Expected Ratios", 
       x = "RepClass", 
       y = "Antisense to Sense Ratio", 
       fill = "Ratio Type") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))



p_values <- ratio_per_gene %>%
  inner_join(observed_ratios, by = "repClass") %>%
  group_by(repClass) %>%
  #Allow p value to be super small
  summarize(p_value = wilcox.test(gene_ratio, observed_ratio, paired = F)$p.value)


plot_data <- ratio_per_gene %>%
  dplyr::select(repClass, gene_ratio) %>%
  mutate(Ratio_Type = "Expected") %>%
  bind_rows(mutate(observed_ratios, 
                   gene_ratio = observed_ratio, 
                   Ratio_Type = "Observed"))

ggboxplot(plot_data, x = "repClass", y = "gene_ratio",
          color = "Ratio_Type", palette = c("#00AFBB", "#E7B800"),
          add = "jitter", shape = "Ratio_Type",
          title = "Comparison of Observed vs Expected Ratios",
          ylab = "Antisense to Sense Ratio", xlab = "RepClass") +
  stat_compare_means(aes(group = Ratio_Type), 
                     method = "wilcox.test", paired = F, 
                     label = "p.signif", label.y = max(plot_data$gene_ratio))


# Create the density plot for each repClass
p <- ggplot(ratio_per_gene, aes(x = gene_ratio)) +
  geom_density(fill = "#00AFBB", alpha = 0.5) +
  geom_vline(data = observed_ratios, aes(xintercept = observed_ratio, color = repClass), 
             linetype = "dashed", size = 1) +
  facet_wrap(~ repClass, scales = "free") +
  labs(title = "Distribution of Expected Ratios with Observed Ratios",
       y = "Density", x = "Antisense to Sense Ratio") +
  theme_minimal() +
  scale_color_manual(values = c("red", "blue", "green", "purple"))

print(p)

p <- ggplot(ratio_per_gene, aes(x = repClass, y = gene_ratio)) +
  geom_violin(fill = "#00AFBB", alpha = 0.5) +
  geom_point(data = observed_ratios, aes(x = repClass, y = observed_ratio), color = "red", size = 3) +
  labs(title = "Boxplot of Expected Ratios with Observed Ratios",
       y = "Antisense to Sense Ratio", x = "RepClass") +
  theme_minimal() +
  scale_y_log10() +
  #Add empirical p-value from combined_ratios
  geom_text(data = combined_ratios, aes(x = repClass, y = 40, label = paste0("p = ", round(p_value, 5))), 
            color = "red", size = 5, vjust = -1) +
  #Resize fonts
  theme(axis.text.x = element_text(color = "grey20", size = 12, angle = 0, hjust = .5, vjust = .5, face = "plain"),
        axis.text.y = element_text(color = "grey20", size = 14, angle = 0, hjust = 1, vjust = 0, face = "plain"),  
        axis.title.x = element_text(color = "grey20", size = 14, angle = 0, hjust = .5, vjust = 0, face = "plain"),
        axis.title.y = element_text(color = "grey20", size = 14, angle = 90, hjust = .5, vjust = .5, face = "plain"),
        plot.title = element_text(color = "grey20", size = 16, angle = 0, hjust = .5, vjust = .5, face = "plain"))
pdf(paste0(plots_path, "boxplot_of_expected_ratios_with_observed_ratios.pdf"), width = 8, height = 6)
print(p)
dev.off()

```

## Cassette exon ratio statistics story

```{r}
library(ggplot2)


# Compute the mean ratio for each repClass
expected_ratio <- count_of_tes_oriented %>%
  group_by(repClass, te_orientation) %>%
  summarise(total_tes = sum(number_of_tes)) %>%
  spread(te_orientation, total_tes) %>%
  mutate(mean_ratio = antisense/sense) %>%
  dplyr::select(repClass, mean_ratio)


# Plot a boxplot
p1 <- ggviolin(
  count_of_tes_oriented, 
  x = "repClass", 
  y = "number_of_tes", 
  fill = "te_orientation",
  color = "black",
  add = "boxplot",
  ylab = "Number of TEs", 
  xlab = "RepClass",
  title = "Expected Number of TEs by Orientation"
) + 
theme_minimal() +
geom_text(data = expected_ratio, 
          aes(x = repClass, y = max(count_of_tes_oriented$number_of_tes, na.rm = TRUE), label = sprintf("Ratio: %.2f", mean_ratio)),
          vjust = -1, 
          color = "black", 
          inherit.aes = FALSE) +
scale_y_log10() +
theme_pubr()

print(p1)







# Plot

observed_ratio_splice <- included_cassette_te_events %>%
  group_by(repClass, repOrientation) %>%
  summarise(total_tes = sum(n)) %>%
  spread(repOrientation, total_tes) %>%
  mutate(mean_ratio = antisense/sense) %>%
  dplyr::select(repClass, mean_ratio)

p2 <- ggplot(included_cassette_te_events, aes(x = repClass, y = n, fill = repOrientation)) +
  geom_bar(stat = "identity", position = "dodge") +
  labs(title = "Observed Number of Spliced TEs by Orientation",
       y = "Number of TEs", x = "RepClass") +
  theme_minimal() +
  geom_text(data = observed_ratio_splice, aes(x = repClass, y = max(count_of_tes_oriented$number_of_tes, na.rm = TRUE), label = sprintf("Ratio: %.2f", mean_ratio)),
          color = "black", inherit.aes = FALSE) +
  scale_y_log10()

print(p2)



# Calculate empirical p-values for each repClass
empirical_p_values <- lapply(unique(ratio_per_gene$repClass), function(rc) {
  expected <- ratio_per_gene$gene_ratio[ratio_per_gene$repClass == rc]
  observed <- observed_ratios$observed_ratio[observed_ratios$repClass == rc]
  
  # Proportion of expected values greater than (or less than, depending on your hypothesis) the observed
  p_value <- sum(expected > observed) / length(expected)
  
  return(c(repClass = rc, empirical_p_value = p_value))
})

empirical_p_values <- do.call(rbind, empirical_p_values)


# Custom function to determine the position to place the p-value annotations
annotate_position <- function(data) {
  max_val <- max(data$gene_ratio, na.rm = TRUE)
  return(max_val - (0.05 * max_val)) # Adjust as needed
}

ggplot(ratio_per_gene, aes(x = gene_ratio)) +
  geom_density(aes(fill = repClass), alpha = 0.5) +
  geom_vline(data = observed_ratios, aes(xintercept = observed_ratio, color = repClass),
             linetype = "dashed", size = 1) +
  facet_wrap(~ repClass, scales = "free") +
  labs(title = "Distribution of Expected TE Gene Orientation Ratios with Observed Ratios",
       y = "Density", x = "Antisense to Sense Ratio") +
  theme_minimal() +
  scale_fill_manual(values = c("blue", "green", "purple", "orange")) +
  scale_color_manual(values = c("blue", "green", "purple", "orange")) +
  geom_text(data = as.data.frame(empirical_p_values),
            aes(label = sprintf("empirical p = %.3f", as.numeric(empirical_p_value))),
            x = Inf, y = Inf, vjust = 2, hjust = 1.5, inherit.aes = FALSE)




# Density plot
p3 <-ggplot(ratio_per_gene, aes(x = gene_ratio)) +
  geom_density(aes(fill = repClass), alpha = 0.5) +
  geom_vline(data = observed_ratios, aes(xintercept = observed_ratio, color = repClass),
             linetype = "dashed", size = 1) +
  facet_wrap(~ repClass, scales = "free") +
  labs(title = "Distribution of Expected TE Gene Orientation Ratios with Spliced TE Ratios (depicted as dashed lines)",
       y = "Density", x = "Antisense to Sense Ratio") +
  theme_minimal() +
  scale_fill_manual(values = c("blue", "green", "purple", "orange")) +
  scale_color_manual(values = c("blue", "green", "purple", "orange")) +
  geom_text(data = as.data.frame(empirical_p_values),
            aes(label = sprintf("empirical p = %.3f", as.numeric(empirical_p_value))),
            x = Inf, y = Inf, vjust = 2, hjust = 1.5, inherit.aes = FALSE)


# Boxplot with observed ratios overlaid
p4 <- ggplot(ratio_per_gene, aes(x = repClass, y = gene_ratio)) +
  geom_boxplot(aes(fill = repClass), alpha = 0.5) +
  geom_point(data = observed_ratios, aes(x = repClass, y = observed_ratio), color = "red", size = 3) +
  labs(title = "Boxplot of Expected TE Gene Orientation Ratios with Spliced TE Ratios (depicted as red circle)",
       y = "Antisense to Sense Ratio", x = "RepClass") +
  theme_minimal() +
  scale_fill_manual(values = c("blue", "green", "purple", "orange")) +
  scale_y_log10() +
  geom_text(data = as.data.frame(empirical_p_values),
            aes(x = repClass, y = Inf, label = sprintf("empirical p = %.3f", as.numeric(empirical_p_value))),
            vjust = 2, inherit.aes = FALSE)



pdf(paste0(plots_path, "cassette_exon_ratio_plots.pdf"), width = 8, height = 5)
print(p1)
print(p2)
print(p3)
print(p4)
dev.off()


```


## Investigate L1 ASP
```{r}
total_AF %>% 
  filter(repClass == "LINE") %>%
  ggplot(aes(x = repName, fill = repOrientation)) +
  geom_bar() +
  theme_minimal() +
  ggtitle("LINEs in Alternative First Exons") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1)) +
  xlab("LINE Name") +
  ylab("Number of transcription start sites")


## Percentages of L1 ASP (Total)
total_AF %>% 
  mutate(AntisenseLINEorNot = ifelse(repClass == "LINE"  & repOrientation == "antisense", "Antisense LINE", "Not Antisense LINE")) %>%
  group_by(AntisenseLINEorNot) %>%
  count() %>%
  ungroup() %>%
  mutate(percent = n / sum(n) * 100)



## Percentages of L1 ASP (Just Lines)
total_AF %>% 
  filter(repClass == "LINE") %>%
  group_by(repOrientation) %>%
  count() %>%
  ungroup() %>%
  mutate(percent = n / sum(n))

# Percentages of L1 ASP in L1HS, L1PA2, L1PA3
total_AF %>% 
  filter(repClass == "LINE", 
         repOrientation == "antisense") %>%
  mutate(Young_L1 = ifelse(repName %in% c("L1HS", "L1PA2", "L1PA3"), "Young L1", "Not Young L1")) %>%
  group_by(Young_L1, repOrientation) %>%
  count() %>%
  ungroup() %>%
  mutate(percent = n / sum(n))



```


```{r}
#Make a logo plot fo the Line elements in the first exons of the genes with alternative TSS
AF_event_grange <- makeGRangesFromDataFrame(total_AF, 
                                            seqnames.field = "seqnames", 
                                            start.field = "start", 
                                            end.field = "end", 
                                            strand.field = "strand", 
                                            keep.extra.columns = TRUE) %>% 
                  mutate(AF_ensg = gene)

# Plot AF events by repFamily
total_AF %>% 
  filter(repClass == "LINE") %>%
  ggplot(aes(x = repFamily, fill = repOrientation)) +
  geom_bar() +
  theme_minimal() +
  ggtitle("LINEs in Alternative First Exons") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1)) +
  xlab("LINE Family") +
  ylab("Number of transcription start sites")

# Plot AF events by repName
pdf(paste0(plots_path, "AF_events_by_repName.pdf"), width = 5, height = 5)
total_AF %>% 
  filter(repClass == "LINE") %>%
  # ORDER LINES BY evoluationary age L1HS, L1PA###, L1PB###, L1MB###, L1MC###, L1MD###, L1ME###, L1MF###, L1MG###, L2, L3, L4
  mutate(repName = factor(repName, levels = rev(c("L1HS", "L1PA2", "L1PA3", "L1PA4", "L1PA14",
                                              "L1PB1",
                                              "L1M3", "L1M5", "L1MB8", "L1ME3A", "L1ME4a",
                                              "L2", "L2a", "L2b", "L2c", "L4_B_Mam")))) %>%
  ggplot(aes(x = repName, fill = repOrientation)) +
  geom_bar() +
  theme_minimal() +
  ggtitle("LINEs in Alternative First Exons") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1)) +
  xlab("LINE Name") +
  ylab("Number of exonization events") +
  scale_fill_manual(values = c("sense" = "#606060", "antisense" = "#D8E1E7")) +
  theme(axis.text.x = element_text(color = "grey20", size = 12, angle = 0, hjust = .5, vjust = .5, face = "plain"),
        axis.text.y = element_text(color = "grey20", size = 14, angle = 0, hjust = 1, vjust = 0, face = "plain"),  
        axis.title.x = element_text(color = "grey20", size = 14, angle = 0, hjust = .5, vjust = 0, face = "plain"),
        axis.title.y = element_text(color = "grey20", size = 14, angle = 90, hjust = .5, vjust = .5, face = "plain"),
        plot.title = element_text(color = "grey20", size = 16, angle = 0, hjust = .5, vjust = .5, face = "plain")) +
  coord_flip()+
  theme_bw()
#Rotate xlab 45 degrees
# theme(axis.text.x = element_text(angle = 45, hjust = .5))
  dev.off()


total_AF_simp <-
total_AF %>% 
  filter(repClass == "LINE") %>%
  plyranges::select(Gene, seqnames, start_1, end_1, start_2, end_2, repName, repClass, repFamily, repOrientation, deltaPSI_TvsN.x)


# count LINEs in each orientation
total_AF %>%
  filter(repClass == "LINE") %>%
  count(repOrientation)


# Get gene IDs of antisense LINEs
total_AF %>%
  filter(repClass == "LINE") %>%
  filter(repOrientation == "antisense") %>%
  arrange(desc(deltaPSI_TvsN.x)) %>%
  distinct(Gene)


#Get nucleotides of first 10bp of promoter sequence
AS_AF_LINES <-
AF_event_grange %>%
  filter(repClass == "LINE")

L1_TSS <- mutate(anchor_start(AS_AF_LINES), width = 300) %>% shift_upstream(150)

# ?plyranges::anchor_start
mutate(anchor_start(AF_event_grange), width = 200) %>% shift_upstream(10) %>% filter(Gene == "GRHL3")

#load hg38 genome
library(BSgenome.Hsapiens.UCSC.hg38)

    
#Get nucleotides of first 10bp of promoter sequence
L1_tss_seq <- getSeq(BSgenome.Hsapiens.UCSC.hg38, L1_TSS)

names(L1_tss_seq) <- paste0(L1_TSS$repName, "_", L1_TSS$Gene, "_", L1_TSS$alternative_transcripts)

#export as fasta
L1_tss_seq %>% 
writeXStringSet(paste0(data_path,"L1_tss.fasta"))


# Motif analysis with biostrings
library(Biostrings)

motif <- "TCTGCAGAGGTTACTGCTG"  # Replace with your actual LINE-1 antisense promoter motif

matches <- vmatchPattern(motif, L1_tss_seq, max.mismatch=3, fixed = "subject") # allow for 1 mismatch, for example


# Define the function to extract matches given the DNA string and the IRanges object.
extractMatches <- function(dna_string, iranges) {
  # Initialize a vector to store the extracted sequences
  extracted <- vector("character", length(iranges))
  
  # Loop through the IRanges object
  for(j in seq_along(iranges)) {
    # Get the start and end of the range
    s <- start(iranges[j])
    e <- end(iranges[j])
    
    # Check if the range is within the bounds of the dna_string
    if(s > 0 && e <= nchar(as.character(dna_string))) {
      # If within bounds, extract the sequence
      extracted[j] <- as.character(subseq(dna_string, s, e))
    } else {
      # If out of bounds, assign NA
      extracted[j] <- NA
    }
  }
  return(extracted)
}

# Now let's loop through each match and each sequence in L1_tss_seq and extract the valid sequences
extracted_seqs <- vector("list", length(matches))
for(i in seq_along(matches)) {
  extracted_seqs[[i]] <- extractMatches(L1_tss_seq[[i]], matches[[i]])
}

library(ggseqlogo)

# Combine all matched sequences into one character vector
combined_seqs <- unlist(extracted_seqs)

# Convert to a DNAStringSet if not already
library(Biostrings)
dna_seqs <- DNAStringSet(na.omit(combined_seqs))

library(ggseqlogo)
ggseqlogo(na.omit(combined_seqs))
# motif



```

#OREGANNO TEST
```{r}
#GRANGE HERE
total_AF

total_AF %>% filter(seqnames == "chr14") %>% dplyr::select(Gene)

#TIBBLE HERE
total_AF_simp

library(biomaRt)

listEnsembl()

ensembl <- useEnsembl(biomart = "regulation")

listDatasets(mart=ensembl)

ensembl <- useEnsembl(biomart = "regulation", dataset = "hsapiens_regulatory_feature")

listFilters(ensembl)

listAttributes(ensembl)

#coordinate to search for transcription factor binding sites
# chr17:28,656,064-28,656,825
values_to_check <- list(chromosome_name = "chr17", start = 28651064, end = 28666825)
coordinates_to_check <- list(chromosome_name = "chr1", start = 1000000, end = 1500000)
default_values <- list(chromosome_name = "chr1", start = 1, end = 1000000)

getBM(attributes = c('chromosome_name', 'chromosome_start', 'chromosome_end', 'activity', 'regulatory_stable_id', 'feature_type_name'),
      filters = c('chromosome_name', 'start', 'end'),
      values = default_values, 
      mart = ensembl)




filters_available <- listFilters(ensembl)

# Choose a filter to explore (adjust as needed)
chosen_filter <- filters_available$name[1]


# Specify the chosen filter and a default value
chosen_filter <- "chromosome_name"
default_value <- "chr1"  # Adjust the default value as needed

# Perform the query with the chosen filter and default value
result <- getBM(
  attributes = c('chromosome_name', 'chromosome_start', 'chromosome_end', 'activity', 'regulatory_stable_id', 'feature_type_name'),
  filters = chosen_filter,
  values = list(chosen_filter = default_value),
  mart = ensembl
)






#Load Oreganno
orgeanno_coords <- read_tsv(paste0(data_path, "Regulatory_annotations/oreganoAll_coords.tsv"))
oreganno_details <- read_tsv(paste0(data_path, "Regulatory_annotations/oregano_details_2.txt"))

oreganno_details %>% filter(attribute == "type") %>% distinct(attrVal)

just_tfbs <-
  oreganno_details %>% 
  group_by(`#id`) %>% 
  filter(any(attrVal == "Transcription Factor Binding Site")) %>%
  filter(attribute == "Gene") %>%
  dplyr::select(id = `#id`, attrVal)

orgeanno_coords_tfbs <-
      orgeanno_coords %>% 
      left_join(just_tfbs)


oreganno_grange <-
orgeanno_coords_tfbs %>% 
  makeGRangesFromDataFrame(seqnames.field = "#chrom", 
                                            start.field = "chromStart", 
                                            end.field = "chromEnd", 
                                            strand.field = "strand", 
                                            keep.extra.columns = T) 


#Get overlaps between oreganno and AF events ~200bp upstream of TSS
upstream_AF_events <-
AF_event_grange %>% 
  filter(repClass == "LINE") %>%
  mutate(shift_upstream(.), width = 200) %>%
  shift_upstream(200)


common_oreganno <-
upstream_AF_events %>%
  find_overlaps(oreganno_grange) %>%
  as_tibble() %>%
  dplyr::select(seqnames, start, end, repClass, repOrientation, attrVal) %>%
  distinct() %>%
  count(attrVal, name = "oreganno_count") %>%
  arrange(desc(oreganno_count))


upstream_AF_events %>%
  find_overlaps(oreganno_grange) %>%  
  filter(id == "OREG1167868")
upstream_AF_events %>%
  find_overlaps(oreganno_grange) %>%  
  filter(id == "OREG0110968")


alu_tss <- 
upstream_AF_events %>%
  filter(repClass == "SINE",
         repFamily == "Alu") %>%
  as_tibble()



test <-
upstream_AF_events %>%
  find_overlaps(oreganno_grange) %>%
  filter(id %in% common_oreganno$id) %>%
  as_tibble() %>%
  dplyr::select(Gene, seqnames, start_1, end_1, start_2, end_2, repName, repClass, repFamily, repOrientation, deltaPSI_TvsN.x, id) %>%
  ungroup() %>%
  arrange(id)
  

```

<!-- ## First Exon Logo plot -->
<!-- ```{r} -->
<!-- # Logic -->
<!-- # For each gene, calculate the distance of the alternative TE TSS from the first exon TSS -->
<!-- # Import gencode TSL1 TSS -->
<!-- gencode_first <- import("~/Desktop/All_suppa_AS_events/data/first_exons_gencode_TSL1.gtf") -->
<!-- # Get 5p end for each gencode start -->
<!-- gencode_first_5p <- mutate(anchor_5p(gencode_first), width = 1) %>% mutate(gencode_ensg = gene_id) -->
<!-- # Prepare alternative TSS events -->
<!-- AF_event_grange <- makeGRangesFromDataFrame(total_AF,  -->
<!--                                             seqnames.field = "seqnames",  -->
<!--                                             start.field = "start",  -->
<!--                                             end.field = "end",  -->
<!--                                             strand.field = "strand",  -->
<!--                                             keep.extra.columns = TRUE) %>%  -->
<!--                   mutate(AF_ensg = gene) -->
<!-- # Find the distance of the alternative TE TSS from the first exon TSS -->
<!-- first_exon_joined <- join_nearest(gencode_first_5p, AF_event_grange, distance = T, suffix = c(".gencode", ".isoTE")) -->


<!-- closest_gene_to_each_te <- join_nearest(AF_event_grange, gencode_first_5p, distance = T, suffix = c(".isoTE", ".gencode")) -->




<!-- # For each gene, get the distance of the TSS  -->
<!-- first_exon_joined %>% -->
<!--   as_tibble() %>% -->
<!--   left_join(total_AF, by = c("seqnames", "start", "end", "strand")) %>% -->
<!--   filter(gencode_ensg %in% closest_gene_to_each_te$gencode_ensg) %>% -->
<!--   group_by(gencode_ensg) %>% -->
<!--   ggplot(aes(x = distance, fill = repClass.y)) + -->
<!--   geom_histogram() + -->
<!--   labs(title = "Distance of Alternative TSS from First Exon TSS", x = "Distance", y = "Count") + -->
<!--   theme_bw() -->



<!-- first_exon_joined %>% -->
<!--   as_tibble() %>% -->
<!--   left_join(total_AF, by = c("seqnames", "start", "end", "strand")) %>% -->
<!--   filter(gencode_ensg %in% closest_gene_to_each_te$gencode_ensg) %>% -->
<!--   colnames() -->



<!-- # Get the number of TSS for each gene -->
<!-- gencode_first_5p %>%  -->
<!--   as_tibble() %>% -->
<!--   arrange(seqnames, start, gene_id) %>% -->
<!--   group_by(gene_id, seqnames) %>% -->
<!--   mutate(prev_distance = start - lag(start), -->
<!--          next_distance = lead(start) - start) %>% -->
<!--   filter((is.na(prev_distance) | prev_distance > 50) & -->
<!--          (is.na(next_distance) | next_distance > 50)) %>% -->
<!--   summarise(count_TSS = n()) %>% -->
<!--   ungroup() %>% -->
<!--   ggplot(aes(x = count_TSS)) + -->
<!--   geom_histogram(binwidth = 1) + -->
<!--   labs(title = "Number of TSS per Gene", x = "Number of TSS", y = "Count") + -->
<!--   theme_bw() -->




<!-- gencode_first_5p %>%  -->
<!--   as_tibble() %>% -->
<!--   group_by(gene_id) %>% -->
<!--   dplyr:: count(name = "count_TSS") %>% -->
<!--   ggplot(aes(x = count_TSS)) + -->
<!--   geom_histogram(binwidth = 1) + -->
<!--   labs(title = "Number of TSS per Gene", x = "Number of TSS", y = "Count") + -->
<!--   theme_bw() -->




<!-- # gencode_first %>% -->
<!-- #   as_tibble() %>% -->
<!-- #   ggplot(aes(x = gene_type)) + -->
<!-- #   geom_bar() -->


<!-- first_exon_joined %>% -->
<!--   as_tibble() %>% -->
<!--   filter(gencode_ensg == AF_ensg) %>% -->
<!--   # mutate(both_present = if_else(gencode_ensg == AF_ensg, TRUE, FALSE)) %>% -->
<!--   ggplot(aes(x = both_present, fill = gene_type)) + -->
<!--   geom_bar() -->

<!-- # For each gene, calculate the distance of the alternative TE TSS from the first exon TSS -->
<!-- simplified_distance_and_repclass <- -->
<!-- first_exon_joined %>%  -->
<!--   as_tibble() %>%  -->
<!--   filter(gencode_ensg == AF_ensg) %>% -->
<!--   arrange(gene_name, distance) %>% -->
<!--   dplyr::select(gene_name, distance, repClass, transcript_type) -->

<!-- # Annotate identical starts -->
<!-- simplified_distance_and_repclass %>% -->
<!--   ggplot(aes(x = distance)) + -->
<!--   geom_density(aes(fill = repClass), alpha = 0.5) + -->
<!--   scale_x_continuous(labels = scales::unit_format(unit = "K", scale = 1e-3)) + -->
<!--   theme_bw() -->




<!-- rmsk_grange <-  -->
<!--   rmsk.gr %>%  -->
<!--   filter(repClass %in% c("LTR", "LINE", "SINE", "DNA")) %>% -->
<!--   mutate(repChr = seqnames(.)) %>% -->
<!--   mutate(repStart = start(.)) %>% -->
<!--   mutate(repEnd = end(.)) %>% -->
<!--   mutate(repStrand = strand) -->


<!-- # Within this window for each gene, count the number of TEs in each orientation within a 100kb window -->
<!-- gencode_first_5p %>% join_by() -->


<!-- # Calculate the ratio of antisense to sense TEs -->
<!-- # Compare this ratio to the observed ratio of antisense to sense TEs in the first exon -->
<!-- # Plot the distribution of expected ratios with the observed ratio overlaid -->


<!-- ``` -->

## Last Exon Ratio expected
```{r}

# For each gene, figure out how many last exons are shortening events vs how many are lengthening events
#Import all gencode transcripts
all_gencode_v30 <- import.gff2(paste0(data_path, "/gencode.sorted.gtf"))


  


# If end of last exon is shorter than gencode TLS 1 isoform end, then it is a shortening event

#Make a logo plot fo the Line elements in the first exons of the genes with alternative TSS
AL_event_grange <- makeGRangesFromDataFrame(total_AL, 
                                            seqnames.field = "seqnames", 
                                            start.field = "start", 
                                            end.field = "end", 
                                            strand.field = "strand", 
                                            keep.extra.columns = TRUE) %>% 
                  mutate(AL_ensg = gene)

# Get the end
TSL1_ends <- mutate(anchor_3p(gencode_last))
AL_event_ends <- mutate(anchor_3p(AL_event_grange))


# Join the two by gene
  AL_event_ends %>% 
    as_tibble() %>%
    filter(repClass %in% c("LINE", "SINE")) %>%
    left_join(as_tibble(TSL1_ends), by = c("gene" = "gene_id")) %>%
    mutate(repOrientation = if_else(repStrand == strand.y, "sense", "antisense")) %>%
    # Calculate distance between the two considering strand - will be upstream and + will be downstream
    mutate(distance_from_canonical = if_else(strand.x == "+", start.x - end.y, end.x - start.y)) %>%
    dplyr::select(gene, distance_from_canonical, repClass, repOrientation) %>%
    filter(abs(distance_from_canonical) < 25000) %>%
    ggplot(aes(x = distance_from_canonical, fill = repOrientation)) +
    geom_histogram(alpha = 0.5) +
    theme_bw() +
    #Scale y to be whole numbers
    scale_y_continuous(breaks = scales::pretty_breaks(n = 10)) +
    #Make X axis scale in K for thousands
    scale_x_continuous(labels = scales::unit_format(unit = "Kilobase", scale = 1e-3)) +
    #Make X axis from -25kb to 25kb
    coord_cartesian(xlim = c(-25000, 25000)) +
    facet_wrap(~repClass, ncol = 2) +
    scale_fill_manual(values = c("sense" = "#606060", "antisense" = "#D8E1E7")) +
    labs(title = "Distance from TSL1 end to alternative TE end", x = "Distance (kb)", y = "Count") +
    scale_color_manual("black")

# Calculate the difference when TE is in REF but not in ALT
  
```

## Find how many last exons terminate in a TE vs splice acceptor in a TE
```{r}
AL_event_grange %>%
  as_tibble() %>%
  # TE_5p = acceptor in TE, TE_3p = Transcription termination in TE
  mutate()


pdf(paste0(plots_path, "/TE_overlap_with_last_exon.pdf"), width = 12, height = 4)
AL_event_grange %>%
  as_tibble() %>%
  mutate(Exon_TE_overlap = case_when(TE_5p == TRUE & TE_3p == TRUE ~ "Both",
                                     TE_5p == TRUE & is.na(TE_3p) ~ "TE intersects 5' end of last exon",
                                     is.na(TE_5p) & TE_3p == TRUE ~ "last exon terminates in TE",
                                     is.na(TE_5p) & is.na(TE_3p) ~ "Neither")) %>%

  ggplot(aes(x = Exon_TE_overlap, fill = repOrientation)) +
  geom_bar(position = "dodge") +
  theme_bw() +
  labs(title = "TE overlap with last exon", x = "Overlap", y = "Count", fill = "Orientation") +
  theme(legend.position = "bottom") +
  facet_wrap(~repClass, ncol = 4) +
  scale_fill_manual(values = c("sense" = "#606060", "antisense" = "#D8E1E7")) +
  # Label totals for each
  stat_count(aes(label = after_stat(count)), geom = "text", position = position_dodge(width = 1), hjust = .5) +
  coord_flip() 
dev.off()



```

#RBP enrichment for as events

```{r}
# Load RBP expression data

load(paste0(data_path, "brca_GeneExp_HTSeq_FPKM.rda"))
rbp_list <- read_tsv(paste0(data_path, "RBP_shortList.txt"))




clin.data <- as.data.frame(SummarizedExperiment::colData(data))
clin.data <- clin.data %>% dplyr::filter(shortLetterCode == "TP")
exp.matrix <- SummarizedExperiment::assay(data) #FPKM matrix
exp.matrix.tp <- exp.matrix[ , colnames(exp.matrix) %in% clin.data$barcode ]


# Junctions of interest
top_jxns <- top_each_SE %>% filter(Gene == "RHOA") %>% pull(junction_bak)

results_df <- data.frame()

# For junction, get the tumor ids
for(jxn in top_jxns){
  tumors_with_event <- 
    top_each_SE %>% 
    filter(junction_bak == jxn) %>% 
    pull(Tumor_ids) %>% 
    str_split(pattern = ",") %>% 
    unlist()
  # for each RBP correlate the RBP expression with the PSI
  for(rbp in unique(rbp_list$ensembl_gene_id)){
    rbp_exp <- exp.matrix.tp[rownames(exp.matrix.tp) == rbp, ]
    if(length(rbp_exp) == 0){
      next
    }
    rbp_exp_with_event <- rbp_exp[names(rbp_exp) %in% tumors_with_event]
    rbp_exp_without_event <- rbp_exp[!(names(rbp_exp) %in% tumors_with_event)]
    # Do a wilcox test
    
    result <- broom::tidy(wilcox.test(rbp_exp_with_event, rbp_exp_without_event))
    result$rbp_id <- rbp
    results_df <- rbind(results_df, result)
  }
}


# get top differential RBPs
top_rbps <- results_df %>% 
  # Adjust p values
  mutate(p.adj = p.adjust(p.value, method = "BH")) %>%
  filter(p.adj < 0.05) %>% 
  arrange(p.adj) %>% 
  head(5)

# Convert the RBP ensembl ids to gene names
library(org.Hs.eg.db)
top_rbps$gene_symbol <- mapIds(org.Hs.eg.db, keys = top_rbps$rbp_id, column = "SYMBOL", keytype = "ENSEMBL")

# Make a dataframe of the top RBP expression with and without the event
rbp_exp_df <- tibble()
for(rbp in top_rbps$rbp_id){
  current_rbp <- top_rbps %>% filter(rbp_id == rbp)
  rbp_exp <- exp.matrix.tp[rownames(exp.matrix.tp) == rbp, ]
  rbp_exp_with_event <- rbp_exp[names(rbp_exp) %in% tumors_with_event]
  rbp_exp_without_event <- rbp_exp[!(names(rbp_exp) %in% tumors_with_event)]
  rbp_exp_df <- rbind(rbp_exp_df, data.frame(rbp_symbol = current_rbp$gene_symbol, rbp_id = rbp, expression = c(rbp_exp_with_event, rbp_exp_without_event), event = c(rep("With event", length(rbp_exp_with_event)), rep("Without event", length(rbp_exp_without_event)))))
}

library(ggridges)
# Make a ridge plot of the top RBP expression with and without the event
pdf(paste0(plots_path, "RBP_expression_with_and_without_event_ridge.pdf"), width = 4, height = 8)
rbp_exp_df %>% 
  ggplot(aes(x = log2(expression), y = rbp_symbol, fill = event)) +
  geom_density_ridges(scale = 3, alpha = 0.5) +
  theme_bw() +
  labs(title = "RBP expression with and without event", x = "Log2 scaled FPKM", y = "RBP") +
  # scale_fill_manual(values = c("With event" = "#D8E1E7", "Without event" = "#606060")) +
  scale_fill_manual(values = c("With event" = "red", "Without event" = "green")) +
  theme(legend.position = "bottom")
dev.off()



```

 #RBP enrichment refactored

```{r}
# Load required libraries
library(dplyr)
library(tidyr)
library(broom)
library(org.Hs.eg.db)
library(ggplot2)
library(ggridges)

# Load RBP expression data
load(paste0(data_path, "brca_GeneExp_HTSeq_FPKM.rda"))
rbp_list <- read_tsv(paste0(data_path, "RBP_shortList.txt"))

clin.data <- as.data.frame(SummarizedExperiment::colData(data))
clin.data <- clin.data %>% dplyr::filter(shortLetterCode == "TP")
exp.matrix <- SummarizedExperiment::assay(data) # FPKM matrix
exp.matrix.tp <- exp.matrix[, colnames(exp.matrix) %in% clin.data$barcode]

# Junctions of interest
top_jxns <- total_SE$junction_bak

results_df <- data.frame()

# For each junction, get the tumor ids and perform computations
for (jxn in top_jxns) {
  tumors_with_event <- 
    total_SE %>% 
    filter(junction_bak == jxn) %>% 
    pull(Tumor_ids) %>% 
    str_split(pattern = ",") %>% 
    unlist()
  
  # For each RBP, correlate the RBP expression with the PSI
  for (rbp in unique(rbp_list$ensembl_gene_id)) {
    rbp_exp <- exp.matrix.tp[rownames(exp.matrix.tp) == rbp, ]
    if (length(rbp_exp) == 0) {
      next
    }
    
    rbp_exp_with_event <- rbp_exp[names(rbp_exp) %in% tumors_with_event]
    rbp_exp_without_event <- rbp_exp[!(names(rbp_exp) %in% tumors_with_event)]
    
    if (length(rbp_exp_with_event) < 10 | length(rbp_exp_without_event) < 10) {
      next
    }
    
    # Perform a Wilcoxon test
    result <- broom::tidy(wilcox.test(rbp_exp_with_event, rbp_exp_without_event))
    result$rbp_id <- rbp
    # Add the number of tumors with and without the event
    result$n_with_event <- length(tumors_with_event)
    result$n_without_event <- length(colnames(exp.matrix.tp)) - length(tumors_with_event)
    # Add the junction to the results dataframe
    result$jxn <- jxn
    # Add the gene symbol of the Junction
    result$jxn_gene_symbol <- total_SE %>% filter(junction_bak == jxn) %>% pull(Gene)
    result$te_type <- total_SE %>% filter(junction_bak == jxn) %>% pull(repFamily)
    result$te_orientation <- total_SE %>% filter(junction_bak == jxn) %>% pull(repOrientation)
    result$jxn_gene_symbol_repClass <- paste0(result$jxn_gene_symbol, "_", result$te_type)
    results_df <- bind_rows(results_df, result)
  }
  # Add the gene symbol
  results_df$rbp_gene_symbol <- mapIds(org.Hs.eg.db, keys = results_df$rbp_id, column = "SYMBOL", keytype = "ENSEMBL")
}

# Get top differential RBPs
top_rbps <- 
  results_df %>% 
  as_tibble() %>%
  group_by(jxn) %>%
  # Adjust p-values for multiple testing
  mutate(p.adj = p.adjust(p.value, method = "BH")) %>%
  filter(p.adj < 0.01) %>% 
  arrange(p.adj) 
  # slice(1:5) %>%
  # ungroup()


# Do a test to see if there is a rbp that is overrepresented for a particular te_type
# for each rbp, get the number of junctions it is associated with

rbp_jxn_filtered <-
  results_df %>% 
  group_by(jxn) %>%
  # Adjust p-values for multiple testing
  mutate(p.adj = p.adjust(p.value, method = "BH")) %>%
  filter(p.adj < 0.001)%>%
  mutate(type_orientation = paste0(te_type, "_", te_orientation))

rbp_jxn_count <- 
  results_df %>% 
  group_by(jxn) %>%
  # Adjust p-values for multiple testing
  mutate(p.adj = p.adjust(p.value, method = "BH")) %>%
  filter(p.adj < 0.01) %>% 
  group_by(rbp_gene_symbol, te_type, te_orientation) %>% 
  summarise(n_jxns = n_distinct(jxn),
            mean_pvalue = mean(p.adj)) %>%
  arrange(desc(n_jxns)) 



# Do a test to see if there is a rbp that is overrepresented for a particular te_type and orientation

just_alu <-
rbp_jxn_filtered %>% filter(te_type == "Alu")

contingency_table <- table(just_alu$rbp_gene_symbol, just_alu$type_orientation)

library(rstatix)

# Initialize empty lists to store results
antisense_results <- list()
sense_results <- list()

# Iterate over each RBP
for (rbp in rownames(contingency_table)) {
    # Extract observed counts for the current RBP in both categories
    observed_counts_antisense <- contingency_table[rbp, "Alu_antisense"]
    observed_counts_sense <- contingency_table[rbp, "Alu_sense"]
    
    # Check if there are enough observed counts to perform the Chi-squared test
    if (length(observed_counts_antisense) >= 2) {
        # Perform Chi-squared test for the current RBP in the "Alu_antisense" category
        chi_squared_result_antisense <- chisq.test(observed_counts_antisense)
        antisense_results[[rbp]] <- chi_squared_result_antisense
    } else {
        antisense_results[[rbp]] <- "Insufficient data for Chi-squared test"
    }
    
    if (length(observed_counts_sense) >= 2) {
        # Perform Chi-squared test for the current RBP in the "Alu_sense" category
        chi_squared_result_sense <- chisq.test(observed_counts_sense)
        sense_results[[rbp]] <- chi_squared_result_sense
    } else {
        sense_results[[rbp]] <- "Insufficient data for Chi-squared test"
    }
}

# Print results
print("Alu_antisense results:")
print(antisense_results)

print("Alu_sense results:")
print(sense_results)



# 
# # Get the top 5 RBPs for each TE type
# top_rbps_te <- 
#   rbp_jxn_count %>% 
#   group_by(te_type) %>% 
#   slice(1:5) %>% 
#   ungroup() %>% 
#   arrange(te_type, desc(n_jxns)) %>%
#   filter(n_jxns > 3)
# 
# 
# rbp_jxn_count %>%
#   filter(rbp_gene_symbol == "HNRNPC")

# Make a dataframe of the top RBP expression with and without the event
rbp_exp_df <- tibble()
# For each junction and then for each RBP, get the expression with and without the event

for (jxn_current in unique(top_rbps$jxn)) {
  current_jxn <- top_rbps %>% filter(jxn == jxn_current)
    for (rbp in current_jxn$rbp_id) {
      current_rbp <- current_jxn %>% filter(rbp_id == rbp)
      # If the RBP is not in the expression matrix, skip it
      if (length(exp.matrix.tp[rownames(exp.matrix.tp) == rbp, ]) == 0) {
        next
      }
      rbp_exp <- exp.matrix.tp[rownames(exp.matrix.tp) == rbp, ]
      rbp_exp_with_event <- rbp_exp[names(rbp_exp) %in% tumors_with_event]
      rbp_exp_without_event <- rbp_exp[!(names(rbp_exp) %in% tumors_with_event)]
      # Don't forget to add the junction to the dataframe
      rbp_exp_df <- bind_rows(rbp_exp_df, data.frame(rbp_repclass = current_rbp$te_type, rbp_repclass = current_rbp$jxn_gene_symbol_repClass, rbp_gene_symbol = current_rbp$rbp_gene_symbol, junction_gene_symbol = current_rbp$jxn_gene_symbol, junction = current_rbp$jxn, rbp_id = rbp, expression = c(rbp_exp_with_event, rbp_exp_without_event), event = c(rep("With event", length(rbp_exp_with_event)), rep("Without event", length(rbp_exp_without_event)))))
    }
}

# Make a ridge plot of the top RBP expression with and without the event
pdf(paste0(plots_path, "RBP_expression_with_and_without_event_ridge_all_se.pdf"), width = 24, height = 24)
rbp_exp_df %>% 
  mutate(log2_expression = log2(expression + 1)) %>%
  ggplot(aes(x = (log2_expression), y = rbp_gene_symbol, fill = event)) +
  geom_density_ridges(scale = 3, alpha = 0.5) +
  theme_bw() +
  labs(title = "RBP expression with and without event", x = "Log2 scaled FPKM", y = "RBP") +
  # Make x axis whole numbers
  scale_fill_manual(values = c("With event" = "red", "Without event" = "green")) +
  theme(legend.position = "bottom") +
  facet_wrap(~rbp_repclass.1, scales = "free_y")
dev.off()





```
