---
title: "PlotFig4"
output: html_document
date: '2023-12-15'
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

```{r}
afgr_data
```

```{r message=FALSE, warning=FALSE, paged.print=FALSE}
library(arrow,lib.loc = '~/userRLib/')
library(ggupset,lib.loc = '~/userRLib/')
library(cowplot)
library(dplyr)
library(tidyr)
library(PCAtools)
library(boot)
library(edgeR)
library(tibble)
library(stringr)
library(ggplot2)
library(optparse)
library(GenomicFeatures)
```

```{r}
prioritization_df <- read_parquet('/oak/stanford/groups/smontgom/tannerj/UDN_LRS_watershed/Fig4Stuff/Lolipop_raw_data.parquet')

lr_sr_df <- read_delim('/oak/stanford/groups/smontgom/tannerj/UDN_LRS_watershed/SR_LR_comparison/LR_SR_overlap.50.bed',
                       delim = '\t',col_names = c('chr_LR','start_LR','end_LR','SV_LR','SVTYPE_LR',
                                                  'chr_SR','start_SR','end_SR','SV_SR','SVTYPE_SR'))
lr_sr_ins_df <- read_delim('/oak/stanford/groups/smontgom/tannerj/UDN_LRS_watershed/SR_LR_comparison/LR_SR.INS.merged.noGT.SUPP2.ID.exploded.tsv',
                       delim = '\t',col_names = c('SV_LR','SV_SR'))

lr_sr_shared <- lr_sr_df %>% pull(var = 'SV_LR') %>% unique()
lr_sr_ins_shared <- lr_sr_ins_df %>% pull(var = 'SV_LR') %>% unique()
lr_sr_shared_updated <- c(lr_sr_shared,lr_sr_ins_shared) %>% unique()

upset_long_vs_short <- prioritization_df_compare_sr_lr %>% 
  mutate(Watershed_prioritized = ifelse((Watershed_posterior_max<=0.6)|
                                          (is.na(Watershed_posterior_max)),FALSE,TRUE),
                             CADDSV_prioritized = ifelse((CADDSV<10)|(is.na(CADDSV)),FALSE,TRUE),
                             HPO_matched = ifelse(is.na(Rank),FALSE,TRUE),
                             recordID = paste(SubjectID,GeneName,SV,sep=':')) %>% dplyr::filter(Watershed_prioritized==T | CADDSV_prioritized == T) %>%
  dplyr::select(Watershed_prioritized,
                CADDSV_prioritized, 
                unique_to_long_read) %>% as_tibble(rownames = 'recordID') %>%
  gather(Filter, Member,-recordID) %>% 
  filter(Member) %>%
  dplyr::select(-Member) %>% group_by(recordID) %>%
  summarize(Filters = list(Filter)) %>% 
  ggplot(aes(x = Filters)) +
    geom_bar() +
    scale_x_upset()
upset_long_vs_short #%>% ggsave(file="Figure4/Upset_watershed_caddsv_lr_sr.svg", width=5, height=4, dpi=300)
```
```{r}
prioritization_df_compare_sr_lr <- prioritization_df %>% mutate(unique_to_long_read = ifelse(SV %in% lr_sr_shared_updated, F, T))
p_upset <- prioritization_df_compare_sr_lr %>% mutate(Watershed_prioritized = ifelse((Watershed_posterior_max<=0.6)|(is.na(Watershed_posterior_max)),FALSE,TRUE),
                             CADDSV_prioritized = ifelse((CADDSV<10)|(is.na(CADDSV)),FALSE,TRUE),
                             HPO_matched = ifelse(is.na(Rank),FALSE,TRUE),
                             recordID = paste(SubjectID,GeneName,SV,sep=':')) %>% dplyr::filter(HPO_matched == TRUE) %>%
  dplyr::select(Watershed_prioritized,
                CADDSV_prioritized) %>% as_tibble(rownames = 'recordID') %>%
  gather(Filter, Member,-recordID) %>% 
  filter(Member) %>%
  dplyr::select(-Member) %>% group_by(recordID) %>%
  summarize(Filters = list(Filter)) %>%
  ggplot(aes(x = Filters)) +
    geom_bar() +
    scale_x_upset()
p_upset
```

```{r}
upset_long_vs_short
```
```{r}
prioritization_df_compare_sr_lr %>% 
  filter(Watershed_posterior_max>=0.6 & unique_to_long_read==T) %>% 
  group_by(SubjectID) %>% summarise(count=n()) %>% pull(count) %>% median()
prioritization_df_compare_sr_lr %>% 
  filter(Watershed_posterior_max>=0.6 & unique_to_long_read==T) %>% dim()
```
```{r}
prioritization_df_compare_sr_lr %>% 
  filter(Watershed_posterior_max>=0.6 ) %>% dim()
```
```{r}
prioritization_df %>% filter(Watershed_posterior_max>0.6) %>% mutate (prioritized_by_phen2gene = ifelse(is.na(Rank),FALSE, TRUE)) %>% dplyr::select(-one_of(c('Top5_or_not','Gene_of_interest_or_not','Comments','Mouse Gene Symbol/ID','Categories','Rank', 'ID','Score','Status'))) %>% write_delim('SuppTable7_Watershed_SV_prioritized_variants.tsv',delim = '\t',na = 'NA',col_names = TRUE)
```
```{r}
metadata <-read_delim('TableS1_with_MappingColumn - Sheet1 (1).tsv')
```
```{r}
prioritization_df %>% filter(Watershed_posterior_max>0.6) %>% mutate (prioritized_by_phen2gene = ifelse(is.na(Rank),FALSE, TRUE)) %>% dplyr::select(-one_of(c('Top5_or_not','Gene_of_interest_or_not','Comments','Mouse Gene Symbol/ID','Categories','Rank', 'ID','Score','Status'))) %>% left_join(mapping_to_deID,by = 'SubjectID') %>% 
  dplyr::select(Sample_name, everything()) %>% dplyr::select(-SubjectID) %>%  write_delim('SuppTable7_Watershed_SV_prioritized_variants_deIDed_updated_ID.tsv',delim = '\t',na = 'NA',col_names = TRUE)
```
```{r}
prioritization_df %>% filter(Watershed_posterior_max>0.6) %>% mutate (prioritized_by_phen2gene = ifelse(is.na(Rank),FALSE, TRUE)) %>% dplyr::select(-one_of(c('Top5_or_not','Gene_of_interest_or_not','Comments','Mouse Gene Symbol/ID','Categories','Rank', 'ID','Score','Status')))
```

```{r}
mapping_to_deID <- metadata %>% dplyr::select(UDNID,Sample_name) %>% dplyr::rename(SubjectID = UDNID)
```

```{r}
prioritization_df %>% filter(Gene_of_interest_or_not==TRUE)
```


```{r}
mapping_to_deID %>% filter(Sample_name %in% LRS_samples$Sample_name)
```
```{r}
metadata %>% filter(INDIVIDUAL =='UDN102380')
```

```{r}
prioritization_df_compare_sr_lr %>% filter(Watershed_posterior_max >= 0.6) %>% ggplot(aes(fill=unique_to_long_read, x=SVTYPE)) + 
    geom_bar(position="dodge", stat="count")
```
# !!!! run this for Figure 4A
```{r}
prioritization_df <- read_parquet('/oak/stanford/groups/smontgom/tannerj/UDN_LRS_watershed/Fig4Stuff/Lolipop_raw_data.parquet')
lolipop_plot_data <- rbind(
prioritization_df %>% group_by(SubjectID) %>% summarise(gene_sv_count = n(), filter = 'rare SVs +/- 100kb',filter_type='WGS-only'),
prioritization_df %>% filter((CADDSV > 10)) %>% group_by(SubjectID) %>% summarise(gene_sv_count = n(), filter = 'rare SVs +/- 100kb, CADD-SV > 10',filter_type='WGS-only'),
prioritization_df %>% filter((abs(Whole_Blood_outlier_zscore)>3)|(abs(Fibroblast_outlier_zscore)>3)) %>% group_by(SubjectID) %>% summarise(gene_sv_count = n(), filter = 'rare SVs +/- 100kb, eOutlier',filter_type='WGS-RNA'),
prioritization_df %>% filter(((abs(Whole_Blood_outlier_zscore)>3)|(abs(Fibroblast_outlier_zscore)>3))&(CADDSV>10)) %>% group_by(SubjectID) %>% summarise(gene_sv_count = n(), filter = 'rare SVs +/- 100kb, eOutlier, CADD-SV > 10',filter_type='WGS-RNA'),
prioritization_df %>% filter((abs(Watershed_posterior_max)>0.6)&((abs(Whole_Blood_outlier_zscore)>3)|(abs(Fibroblast_outlier_zscore)>3))) %>% group_by(SubjectID) %>% summarise(gene_sv_count = n(), filter = 'rare SVs +/- 100kb, eOutlier, Watershed max > 0.6',filter_type='WGS-RNA'),
prioritization_df %>% filter(!is.na(Rank)) %>% group_by(SubjectID) %>% summarise(gene_sv_count = n(), filter = 'rare SVs +/- 100kb, HPO match',filter_type='WGS-HPO'),
prioritization_df %>% filter((!is.na(Rank))&((abs(Whole_Blood_outlier_zscore)>3)|(abs(Fibroblast_outlier_zscore)>3))) %>% group_by(SubjectID) %>% summarise(gene_sv_count = n(), filter = 'rare SVs +/- 100kb, eOutlier, HPO match',filter_type='WGS-RNA-HPO'),
prioritization_df %>% filter((!is.na(Rank))&(CADDSV>10)) %>% group_by(SubjectID) %>% summarise(gene_sv_count = n(), filter = 'rare SVs +/- 100kb, HPO match, CADD-SV > 10',filter_type='WGS-HPO'),
prioritization_df %>% filter((!is.na(Rank))&(CADDSV>10)&((abs(Whole_Blood_outlier_zscore)>3)|(abs(Fibroblast_outlier_zscore)>3))) %>% group_by(SubjectID) %>% summarise(gene_sv_count = n(), filter = 'rare SVs +/- 100kb, eOutlier, HPO match, CADD-SV > 10',filter_type='WGS-RNA-HPO'),
prioritization_df %>% filter((!is.na(Rank))&(Watershed_posterior_max>0.6)&((abs(Whole_Blood_outlier_zscore)>3)|(abs(Fibroblast_outlier_zscore)>3))) %>% group_by(SubjectID) %>% summarise(gene_sv_count = n(), filter = 'rare SVs +/- 100kb, eOutlier, HPO match, Watershed max > 0.6',filter_type='WGS-RNA-HPO'))
orders <- rev(c('rare SVs +/- 100kb','rare SVs +/- 100kb, CADD-SV > 10','rare SVs +/- 100kb, HPO match','rare SVs +/- 100kb, HPO match, CADD-SV > 10','rare SVs +/- 100kb, eOutlier','rare SVs +/- 100kb, eOutlier, CADD-SV > 10', 'rare SVs +/- 100kb, eOutlier, Watershed max > 0.6', 'rare SVs +/- 100kb, eOutlier, HPO match','rare SVs +/- 100kb, eOutlier, HPO match, CADD-SV > 10', 'rare SVs +/- 100kb, eOutlier, HPO match, Watershed max > 0.6' ))

e <- ggplot(lolipop_plot_data, aes(x = factor(filter,levels = orders), y = log2(gene_sv_count), color=factor(filter_type)))
# Stripcharts with summary statistics
# Change color by dose groups
e1 <- e + geom_jitter(aes(color = filter_type), position = position_jitter(0.2), size = 1.2) + stat_summary(fun.data="mean_sdl",  fun.args = list(mult=1),color = "red",alpha=0.5) + coord_flip() + theme(aspect.ratio=1/1) + theme(legend.position="bottom") + ggtitle('Number of SV-gene pairs prioritized') + scale_color_manual(values = custom_colors) + labs(x='Rare disease variant filters') +  scale_y_continuous(breaks = c(0, 1, 2, 4, 8), labels = c( "1", "2", "4", "16", "256")) + 
  expand_limits(y=0) + 
  labs(color = 'group') + ylab('gene-variant pairs') + 
  guides(colour=guide_legend(nrow=2))

e2 <- ggplot(lolipop_plot_data, aes(y = factor(filter,levels = orders))) + geom_bar() + labs(x = 'Patient Count') + theme(aspect.ratio=3/1,
                                                                                                                          axis.text.y=element_blank(),
                                                                                                                          axis.title.y=element_blank(),
                                                                                                                          axis.ticks.y=element_blank(),
                                                                                                                          plot.margin=margin(l=-1, unit="cm")) 
figure5a <- ggarrange(e1,e2,nrow = 1,widths = c(10,3),heights = 10)
figure5a %>% ggsave(file="Figure4/filter_plot.svg", width=11, height=5, dpi=300)

```

```{r}
p_upset <- prioritization_df %>% mutate(Watershed_prioritized = ifelse((Watershed_posterior_max<=0.5)|(is.na(Watershed_posterior_max)),FALSE,TRUE),
                             CADDSV_prioritized = ifelse((CADDSV<10)|(is.na(CADDSV)),FALSE,TRUE),
                             HPO_matched = ifelse(is.na(Rank),FALSE,TRUE),
                             recordID = paste(SubjectID,GeneName,SV,sep=':')) %>% dplyr::filter(HPO_matched == TRUE) %>%
  dplyr::select(Watershed_prioritized,
                CADDSV_prioritized) %>% as_tibble(rownames = 'recordID') %>%
  gather(Filter, Member,-recordID) %>% 
  filter(Member) %>%
  dplyr::select(-Member) %>% group_by(recordID) %>%
  summarize(Filters = list(Filter)) %>%
  ggplot(aes(x = Filters)) +
    geom_bar() +
    scale_x_upset()
p_upset %>% ggsave(file="Upset_watershed_caddsv.png", width=5, height=4, dpi=300)
```
```{r}
percentiles <- seq(0.5, 0.99, by = 0.01)
watershed_percentiles <- prioritization_df %>% summarise(across(Watershed_posterior_max, ~quantile(., probs = percentiles)))
caddsv_percentiles <- prioritization_df %>% filter(!is.na(CADDSV)) %>% summarise(across(CADDSV, ~quantile(., probs = percentiles)))
```
```{r}
intersection_col <- c()
intersection_hpo_match_col <- c()
union_col <- c()
union_hpo_matched_col <- c()
for (percentile in percentiles) {
    # Perform some operation with each value
    watershed_percentile <- prioritization_df %>% pull(Watershed_posterior_max) %>% quantile(., probs = percentile)
    caddsv_percentile <- prioritization_df %>% filter(!is.na(CADDSV)) %>% pull(CADDSV) %>% quantile(., probs = percentile)
    intersection <-prioritization_df %>% filter(Watershed_posterior_max>watershed_percentile &
                                 (CADDSV>caddsv_percentile&!is.na(CADDSV))) %>% nrow()
    union <- prioritization_df %>% filter(Watershed_posterior_max>watershed_percentile |
                                 (CADDSV>caddsv_percentile&!is.na(CADDSV))) %>% nrow()
    intersection_hpo_matched <-prioritization_df %>% filter(Watershed_posterior_max>watershed_percentile &
                                 (CADDSV>caddsv_percentile&!is.na(CADDSV)) & !is.na(Rank)) %>% nrow()
    union_hpo_matched <-prioritization_df %>% filter(Watershed_posterior_max>watershed_percentile |
                                 (CADDSV>caddsv_percentile&!is.na(CADDSV)) & !is.na(Rank)) %>% nrow()
    intersection_col <- c(intersection_col,intersection)
    intersection_hpo_match_col <- c(intersection_hpo_match_col,intersection_hpo_matched)
    union_col <- c(union_col,union)
    union_hpo_matched_col <- c(union_hpo_matched_col,union_hpo_matched)
}
overlap_by_percentile <- data_frame(percentile = percentiles,
                                    intersection = intersection_col,
                                    intersection_hpo_match = intersection_hpo_match_col,
                                    union = union_col, 
                                    union_hpo_matched = union_hpo_matched_col)
```
```{r}
overlap_by_percentile %>% mutate(jaccard = intersection/union, jaccard_hpo_matched = intersection_hpo_match/union_hpo_matched) %>% ggplot(aes(x = percentile,y = jaccard)) + geom_point()
```

```{r}
prioritization_df %>% mutate(Watershed_prioritized = ifelse((Watershed_posterior_max<=0.5)|(is.na(Watershed_posterior_max)),FALSE,TRUE),
                             CADDSV_prioritized = ifelse((CADDSV<10)|(is.na(CADDSV)),FALSE,TRUE),
                             HPO_matched = ifelse(is.na(Rank),FALSE,TRUE),
                             recordID = paste(SubjectID,GeneName,SV,sep=':')) %>% dplyr::filter(HPO_matched == F) %>% 
```

```{r}
caddsv_percentiles
```

```{r}
tibble(percentile = percentiles, watershed_threshold = watershed_percentiles %>% pull(Watershed_posterior_max),
       caddsv_threshold = caddsv_percentiles %>% pull(CADDSV))
```

```{r}
prioritization_df_UDN258063 <- prioritization_df %>% filter(SubjectID=='UDN258063')
prioritization_df_UDN886476 <- prioritization_df %>% filter(SubjectID=='UDN886476')

```
```{r}
prioritization_df_UDN886476
```

```{r}
prioritization_df_UDN258063 %>% filter(GeneName=='ENSG00000151327')
```

```{r}
lolipop_plot_data <- rbind(
prioritization_df_UDN258063 %>% group_by(SubjectID) %>% summarise(gene_sv_count = n(),
                                                                  contain_diagnostic_variant = if_else(('617_0_Sniffles2.DEL.2FC7SD' %in% SV)&('617_0_Sniffles2.DEL.2FD3SD' %in% SV),TRUE,FALSE),
                                                                  filter = 'rare SVs +/- 100kb', 
                                                                  filter_type = 'WGS-only'),
prioritization_df_UDN258063 %>% filter((CADDSV > 10)) %>% group_by(SubjectID) %>% summarise(gene_sv_count = n(),
                                                                                            contain_diagnostic_variant = if_else(('617_0_Sniffles2.DEL.2FC7SD' %in% SV)&('617_0_Sniffles2.DEL.2FD3SD' %in% SV),TRUE,FALSE), 
                                                                                            filter = 'rare SVs +/- 100kb, CADD-SV > 10', 
                                                                                            filter_type = 'WGS-only'),
prioritization_df_UDN258063 %>% filter((abs(Whole_Blood_outlier_zscore)>3)|(abs(Fibroblast_outlier_zscore)>3)) %>% group_by(SubjectID) %>% summarise(gene_sv_count = n(),
            contain_diagnostic_variant = if_else(('617_0_Sniffles2.DEL.2FC7SD' %in% SV)&('617_0_Sniffles2.DEL.2FD3SD' %in% SV),TRUE,FALSE), filter = 'rare SVs +/- 100kb, eOutlier', filter_type = 'WGS-RNA'),
prioritization_df_UDN258063 %>% filter(((abs(Whole_Blood_outlier_zscore)>3)|(abs(Fibroblast_outlier_zscore)>3))&(CADDSV>10)) %>% group_by(SubjectID) %>%
  summarise(gene_sv_count = n(),
            contain_diagnostic_variant = if_else(('617_0_Sniffles2.DEL.2FC7SD' %in% SV)&('617_0_Sniffles2.DEL.2FD3SD' %in% SV),TRUE,FALSE), 
            filter = 'rare SVs +/- 100kb, eOutlier, CADD-SV > 10', filter_type = 'WGS-RNA'),
prioritization_df_UDN258063 %>% filter((abs(Watershed_posterior_max)>0.6)) %>% group_by(SubjectID) %>% summarise(gene_sv_count = n(),
            contain_diagnostic_variant = if_else(('617_0_Sniffles2.DEL.2FC7SD' %in% SV)&('617_0_Sniffles2.DEL.2FD3SD' %in% SV),TRUE,FALSE), filter = 'rare SVs +/- 100kb, Watershed max > 0.6', filter_type = 'WGS-RNA'),
prioritization_df_UDN258063 %>% filter(!is.na(Rank)) %>% group_by(SubjectID) %>% summarise(gene_sv_count = n(),
            contain_diagnostic_variant = if_else(('617_0_Sniffles2.DEL.2FC7SD' %in% SV)&('617_0_Sniffles2.DEL.2FD3SD' %in% SV),TRUE,FALSE), 
            filter = 'rare SVs +/- 100kb, HPO match',filter_type='WGS-HPO'),
prioritization_df_UDN258063 %>% filter((!is.na(Rank))&((abs(Whole_Blood_outlier_zscore)>3)|(abs(Fibroblast_outlier_zscore)>3))) %>% group_by(SubjectID) %>% summarise(gene_sv_count = n(),
            contain_diagnostic_variant = if_else(('617_0_Sniffles2.DEL.2FC7SD' %in% SV)&('617_0_Sniffles2.DEL.2FD3SD' %in% SV),TRUE,FALSE), filter = 'rare SVs +/- 100kb, eOutlier, HPO match',filter_type='WGS-RNA-HPO'),
prioritization_df_UDN258063 %>% filter((!is.na(Rank))&(CADDSV>10)) %>% group_by(SubjectID) %>% summarise(gene_sv_count = n(),
            contain_diagnostic_variant = if_else(('617_0_Sniffles2.DEL.2FC7SD' %in% SV)&('617_0_Sniffles2.DEL.2FD3SD' %in% SV),TRUE,FALSE), filter = 'rare SVs +/- 100kb, HPO match, CADD-SV > 10',filter_type='WGS-HPO'),
prioritization_df_UDN258063 %>% filter((!is.na(Rank))&(CADDSV>10)&((abs(Whole_Blood_outlier_zscore)>3)|(abs(Fibroblast_outlier_zscore)>3))) %>% group_by(SubjectID) %>% summarise(gene_sv_count = n(),
            contain_diagnostic_variant = if_else(('617_0_Sniffles2.DEL.2FC7SD' %in% SV)&('617_0_Sniffles2.DEL.2FD3SD' %in% SV),TRUE,FALSE), filter = 'rare SVs +/- 100kb, eOutlier, HPO match, CADD-SV > 10',filter_type='WGS-RNA-HPO'),
prioritization_df_UDN258063 %>% filter((!is.na(Rank))&(Watershed_posterior_max>0.6)) %>% group_by(SubjectID) %>% summarise(gene_sv_count = n(),
            contain_diagnostic_variant = if_else(('617_0_Sniffles2.DEL.2FC7SD' %in% SV)&('617_0_Sniffles2.DEL.2FD3SD' %in% SV),TRUE,FALSE), filter = 'rare SVs +/- 100kb, HPO match, Watershed max > 0.6',filter_type='WGS-RNA-HPO'))
orders <- rev(c('rare SVs +/- 100kb','rare SVs +/- 100kb, CADD-SV > 10','rare SVs +/- 100kb, HPO match','rare SVs +/- 100kb, HPO match, CADD-SV > 10','rare SVs +/- 100kb, eOutlier','rare SVs +/- 100kb, eOutlier, CADD-SV > 10', 'rare SVs +/- 100kb, eOutlier, Watershed max > 0.6', 'rare SVs +/- 100kb, eOutlier, HPO match','rare SVs +/- 100kb, eOutlier, HPO match, CADD-SV > 10', 'rare SVs +/- 100kb, eOutlier, HPO match, Watershed max > 0.6' ))
```

```{r}
fig4_udn258063 <- lolipop_plot_data %>% 
ggplot( aes(x=factor(filter), y=gene_sv_count, color=as.factor(filter_type)), hide) +
  geom_segment( aes(x=factor(filter), xend=factor(filter), y=0, yend=gene_sv_count)) +
  geom_point( aes(x=factor(filter),y=gene_sv_count,shape=contain_diagnostic_variant), size=5) + 
  scale_y_continuous(trans = "log", breaks = c(1, 10, 100), labels = c( "1", "10", "100")) + 
  theme_light() +
  coord_flip() +
  scale_color_manual(values = custom_colors) +
  theme(
    panel.grid.major.y = element_blank(),
    panel.border = element_blank(),
    axis.ticks.y = element_blank()
  ) + labs(shape = "Contain Diagnostic Variant", color = "Filter type") + 
  labs( y = "Gene-SV pair count", x = "") 
fig4_udn258063 %>% ggsave(file="Figure4/udn258063.svg", width=10, height=5, dpi=300)

```
```{r}
lolipop_plot_data_vntr <- rbind(
prioritization_df_UDN886476 %>% group_by(SubjectID) %>% summarise(gene_sv_count = n(),
                                                                  contain_diagnostic_variant = if_else(('chr5:177554489:CGC' %in% SV | '599_0_Sniffles2.INS.1325S4' %in% SV),TRUE,FALSE),
                                                                  filter = 'rare SVs +/- 100kb', 
                                                                  filter_type = 'WGS-only'),
prioritization_df_UDN886476 %>% filter((CADDSV > 10)) %>% group_by(SubjectID) %>% summarise(gene_sv_count = n(),
                                                                                            contain_diagnostic_variant = if_else(('chr5:177554489:CGC' %in% SV | '599_0_Sniffles2.INS.1325S4' %in% SV),TRUE,FALSE), 
                                                                                            filter = 'rare SVs +/- 100kb, CADD-SV > 10', 
                                                                                            filter_type = 'WGS-only'),
prioritization_df_UDN886476 %>% filter((abs(Whole_Blood_outlier_zscore)>3)|(abs(Fibroblast_outlier_zscore)>3)) %>% group_by(SubjectID) %>% summarise(gene_sv_count = n(),
            contain_diagnostic_variant = if_else(('chr5:177554489:CGC' %in% SV | '599_0_Sniffles2.INS.1325S4' %in% SV),TRUE,FALSE), filter = 'rare SVs +/- 100kb, eOutlier', filter_type = 'WGS-RNA'),
prioritization_df_UDN886476 %>% filter(((abs(Whole_Blood_outlier_zscore)>3)|(abs(Fibroblast_outlier_zscore)>3))&(CADDSV>10)) %>% group_by(SubjectID) %>%
  summarise(gene_sv_count = n(),
            contain_diagnostic_variant = if_else(('chr5:177554489:CGC' %in% SV | '599_0_Sniffles2.INS.1325S4' %in% SV),TRUE,FALSE), 
            filter = 'rare SVs +/- 100kb, eOutlier, CADD-SV > 10', filter_type = 'WGS-RNA'),
prioritization_df_UDN886476 %>% filter((abs(Watershed_posterior_max)>0.6)) %>% group_by(SubjectID) %>% summarise(gene_sv_count = n(),
            contain_diagnostic_variant = if_else(('chr5:177554489:CGC' %in% SV | '599_0_Sniffles2.INS.1325S4' %in% SV),TRUE,FALSE), filter = 'rare SVs +/- 100kb, Watershed max > 0.6', filter_type = 'WGS-RNA'),
prioritization_df_UDN886476 %>% filter(!is.na(Rank)) %>% group_by(SubjectID) %>% summarise(gene_sv_count = n(),
            contain_diagnostic_variant = if_else(('chr5:177554489:CGC' %in% SV | '599_0_Sniffles2.INS.1325S4' %in% SV),TRUE,FALSE), 
            filter = 'rare SVs +/- 100kb, HPO match',filter_type='WGS-HPO'),
prioritization_df_UDN886476 %>% filter((!is.na(Rank))&((abs(Whole_Blood_outlier_zscore)>3)|(abs(Fibroblast_outlier_zscore)>3))) %>% group_by(SubjectID) %>% summarise(gene_sv_count = n(),
            contain_diagnostic_variant = if_else(('chr5:177554489:CGC' %in% SV | '599_0_Sniffles2.INS.1325S4' %in% SV),TRUE,FALSE), filter = 'rare SVs +/- 100kb, eOutlier, HPO match',filter_type='WGS-RNA-HPO'),
prioritization_df_UDN886476 %>% filter((!is.na(Rank))&(CADDSV>10)) %>% group_by(SubjectID) %>% summarise(gene_sv_count = n(),
            contain_diagnostic_variant = if_else(('chr5:177554489:CGC' %in% SV | '599_0_Sniffles2.INS.1325S4' %in% SV),TRUE,FALSE), filter = 'rare SVs +/- 100kb, HPO match, CADD-SV > 10',filter_type='WGS-HPO'),
prioritization_df_UDN886476 %>% filter((!is.na(Rank))&(CADDSV>10)&((abs(Whole_Blood_outlier_zscore)>3)|(abs(Fibroblast_outlier_zscore)>3))) %>% group_by(SubjectID) %>% summarise(gene_sv_count = n(),
            contain_diagnostic_variant = if_else(('chr5:177554489:CGC' %in% SV | '599_0_Sniffles2.INS.1325S4' %in% SV),TRUE,FALSE), filter = 'rare SVs +/- 100kb, eOutlier, HPO match, CADD-SV > 10',filter_type='WGS-RNA-HPO'),
prioritization_df_UDN886476 %>% filter((!is.na(Rank))&(Watershed_posterior_max>0.6)) %>% group_by(SubjectID) %>% summarise(gene_sv_count = n(),
            contain_diagnostic_variant = if_else(('chr5:177554489:CGC' %in% SV | '599_0_Sniffles2.INS.1325S4' %in% SV),TRUE,FALSE), filter = 'rare SVs +/- 100kb, HPO match, Watershed max > 0.6',filter_type='WGS-RNA-HPO'))
orders <- rev(c('rare SVs +/- 100kb','rare SVs +/- 100kb, CADD-SV > 10','rare SVs +/- 100kb, HPO match','rare SVs +/- 100kb, HPO match, CADD-SV > 10','rare SVs +/- 100kb, eOutlier','rare SVs +/- 100kb, eOutlier, CADD-SV > 10', 'rare SVs +/- 100kb, eOutlier, Watershed max > 0.6', 'rare SVs +/- 100kb, eOutlier, HPO match','rare SVs +/- 100kb, eOutlier, HPO match, CADD-SV > 10', 'rare SVs +/- 100kb, eOutlier, HPO match, Watershed max > 0.6' ))
```

```{r}
custom_colors <- c("WGS-only" = "#FFBC7B", "WGS-RNA" = "#1E88E5", "WGS-HPO" = "#004D40", "WGS-RNA-HPO" = "#D81B60")

fig4_udn886476 <- lolipop_plot_data_vntr %>% 
ggplot( aes(x=factor(filter), y=gene_sv_count, color=as.factor(filter_type)), hide) +
  geom_segment( aes(x=factor(filter), xend=factor(filter), y=0, yend=gene_sv_count)) +
  geom_point( aes(x=factor(filter),y=gene_sv_count,shape=contain_diagnostic_variant), size=5) + 
  scale_y_continuous(trans = "log", breaks = c(1, 10, 100), labels = c( "1", "10", "100")) + 
  theme_light() +
  coord_flip() +
  theme(
    panel.grid.major.y = element_blank(),
    panel.border = element_blank(),
    axis.ticks.y = element_blank()
  ) + scale_color_manual(values = custom_colors) +
  labs(shape = "Contain Diagnostic Variant", color = "Filter type") + 
  labs( y = "Gene-SV pair count", x = "")
fig4_udn886476 %>% ggsave(file="Figure4/udn886476.svg", width=10, height=5, dpi=300)
```
```{r}
prioritization_df_UDN886476 %>% filter(((abs(Whole_Blood_outlier_zscore)>3)|(abs(Fibroblast_outlier_zscore)>3))&(CADDSV>10)) %>% group_by(SubjectID) %>%
  summarise(gene_sv_count = n(),
            contain_diagnostic_variant = if_else(('chr5:177554489:CGC' %in% SV | '599_0_Sniffles2.INS.1325S4' %in% SV),TRUE,FALSE), 
            filter = 'rare SVs +/- 100kb, eOutlier, CADD-SV > 10', filter_type = 'WGS-RNA')
```

```{r}
prioritization_df_UDN886476 %>% filter((abs(Watershed_posterior_max)>0.6)) %>% group_by(SubjectID)
```

```{r}
# does prioritization between cadd and watershed converge when closer to genes. 
prioritization_df %>% filter(!is.na(CADDSV)) %>% ggplot(aes(x=Watershed_posterior_max, y=CADDSV)) + geom_point()
# plot variant distance to genes as x axis, and y axis is differences in prioritization rank between w and c. 
```
```{r}
annotations <- read.delim('/oak/stanford/groups/smontgom/tannerj/UDN_LRS_watershed/Watershed-SV/UDN-protein-lincRNA-100k-ADRC-all-position-aware-2.0/tissue_specific_annotations_pred_no_vep_no_muscle_newimpute_testing_data.tsv',sep = '\t')
```
```{r}
additional_annotations <- annotations %>% separate(GeneName, into = c("GeneName", "SV"), sep = ":", remove = TRUE) %>% dplyr::select(c(SubjectID, GeneName, SV, exon_variant, nc_dist2gene))
```
```{r}
revision_data
```

```{r}
revision_data <- prioritization_df %>% 
  filter(!is.na(CADDSV)) %>% 
  mutate(
    Watershed_rank = rank(Watershed_posterior_max, ties.method = 'first', ),
    CADDSV_rank = rank(CADDSV, ties.method = 'first'),
    Watershed_prioritized = ifelse(Watershed_posterior_max>0.6,T,F),
    CADDSV_prioritized = ifelse(CADDSV>10,T,F)) %>% 
  mutate(
    delta_rank = Watershed_rank-CADDSV_rank) %>% 
  left_join(additional_annotations, by = c('SubjectID', 'GeneName', 'SV')) %>% mutate(
    distance_categories = case_when(nc_dist2gene == 0 & exon_variant == 1 ~ 'coding variant',
                                    nc_dist2gene == 0 & exon_variant == 0 ~ 'intron variant',
                                    nc_dist2gene > 0 & nc_dist2gene < 1000 ~ '1-1kb variant',
                                    nc_dist2gene > 1000 & nc_dist2gene < 10000 ~ '1kb-10kb variant',
                                    nc_dist2gene > 10000 & nc_dist2gene < 100000 ~ '10-100kb variant'))
```
```{r}
revision_data
```
```{r}
df_counts
```

```{r}
df_counts <- revision_data %>%
  group_by(distance_categories) %>%
  summarize(count = n())
revision_data %>% ggplot( aes(x = fct_relevel(distance_categories,
                                              "gene body variant", 
                                              "1-1kb variant",
                                              "1kb-10kb variant",
                                              "10-100kb variant"), y = delta_rank)) +
  geom_violin() +
  geom_boxplot(width = 0.1, fill = "white", color = "black", alpha = 0.7) +
  stat_summary(fun = function(x) min(x) - 0.2 * diff(range(x)), 
               geom = "text", 
               aes(label = paste("n =", df_counts$count[df_counts$distance_categories == unique(x)])),
               vjust = 1.5) +
  labs(x = "Exon Variant", y = "Delta Rank", 
       title = "Distribution of Watershed Rank - CADD-SV Rank by Variant Distance Category") +
  theme_minimal()
```

```{r}
revision_data %>% ggplot( aes(x = factor(exon_variant), y = delta_rank)) +
  geom_violin() +
  geom_boxplot(width = 0.1, fill = "white", color = "black", alpha = 0.7) +
  labs(x = "Exon Variant", y = "Delta Rank", 
       title = "Distribution of Delta Rank by Exon Variant") +
  theme_minimal()
```
```{r}
revision_data_w_quantile
```

```{r}
bin_thresholds <- c(0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9)
revision_data_w_quantile <- revision_data %>% mutate(Watershed_quantile = case_when(
      Watershed_posterior_max <= quantile(Watershed_posterior_max, bin_thresholds[1]) ~ "0-10%",
      Watershed_posterior_max <= quantile(Watershed_posterior_max, bin_thresholds[2]) ~ "10-20%",
      Watershed_posterior_max <= quantile(Watershed_posterior_max, bin_thresholds[3]) ~ "20-30%",
      Watershed_posterior_max <= quantile(Watershed_posterior_max, bin_thresholds[4]) ~ "30-40%",
      Watershed_posterior_max <= quantile(Watershed_posterior_max, bin_thresholds[5]) ~ "40-50%",
      Watershed_posterior_max <= quantile(Watershed_posterior_max, bin_thresholds[6]) ~ "50-60%",
      Watershed_posterior_max <= quantile(Watershed_posterior_max, bin_thresholds[7]) ~ "60-70%",
      Watershed_posterior_max <= quantile(Watershed_posterior_max, bin_thresholds[8]) ~ "70-80%",
      Watershed_posterior_max <= quantile(Watershed_posterior_max, bin_thresholds[9]) ~ "80-90%",
      TRUE ~ "90-100%"),
                                                     CADDSV_quantile = case_when(
      CADDSV <= quantile(CADDSV, bin_thresholds[1]) ~ "0-10%",
      CADDSV <= quantile(CADDSV, bin_thresholds[2]) ~ "10-20%",
      CADDSV <= quantile(CADDSV, bin_thresholds[3]) ~ "20-30%",
      CADDSV <= quantile(CADDSV, bin_thresholds[4]) ~ "30-40%",
      CADDSV <= quantile(CADDSV, bin_thresholds[5]) ~ "40-50%",
      CADDSV <= quantile(CADDSV, bin_thresholds[6]) ~ "50-60%",
      CADDSV <= quantile(CADDSV, bin_thresholds[7]) ~ "60-70%",
      CADDSV <= quantile(CADDSV, bin_thresholds[8]) ~ "70-80%",
      CADDSV <= quantile(CADDSV, bin_thresholds[9]) ~ "80-90%",
      TRUE ~ "90-100%"),)
```

```{r}
category_order <- c('coding variant','intron variant','1-1kb variant','1kb-10kb variant', '10-100kb variant')
p_compare <- revision_data_w_quantile %>% dplyr::select(nc_dist2gene, distance_categories, Watershed_prioritized, CADDSV_prioritized) %>% pivot_longer(cols = c(Watershed_prioritized, CADDSV_prioritized), 
               names_to = "Model", 
               values_to = "Prioritized") %>% mutate(Model = str_remove(Model, "_prioritized$")) %>% 
  group_by(Model, distance_categories,Prioritized) %>% summarise(count = n()) %>%
  group_by(Model,distance_categories) %>% mutate(proportion = count / sum(count)) %>%
  ungroup() %>% 
  ggplot(aes(x = Prioritized, y = proportion, fill = factor(distance_categories,levels = category_order))) +
  geom_bar(stat = "identity", position = position_dodge(width = 0.9)) + 
  scale_fill_discrete() +
  facet_grid(rows = vars(Model)) +
  labs(x = "Prioritized or not", y = "fraction", fill = 'distance categories') +
  theme_minimal() +
  ggtitle("Grouped Bar Plot")
ggsave('Watershed_vs_CADD_distance.prioritized_or_not.png',p_compare,height = 8,width = 16)
```
```{r}
category_order <- c('coding variant','intron variant','1-1kb variant','1kb-10kb variant', '10-100kb variant')
p_compare <- revision_data_w_quantile %>% dplyr::select(nc_dist2gene, distance_categories, Watershed_prioritized, CADDSV_prioritized) %>% pivot_longer(cols = c(Watershed_prioritized, CADDSV_prioritized), 
               names_to = "Model", 
               values_to = "Prioritized") %>% mutate(Model = str_remove(Model, "_prioritized$")) %>% 
  group_by(Model, distance_categories,Prioritized) %>% summarise(count = n()) %>%
  group_by(Model,distance_categories) %>% mutate(proportion = count / sum(count)) %>%
  ungroup() %>% filter(Model == 'Watershed' & Prioritized == TRUE) %>% 
  ggplot(aes(x = Prioritized, y = proportion, fill = factor(distance_categories,levels = category_order))) +
  geom_bar(stat = "identity", position = position_dodge(width = 0.9)) + 
  scale_fill_discrete() +
  labs(x = "Prioritized or not", y = "fraction", fill = 'distance categories') +
  theme_minimal() +
  ggtitle("Proportion of prioritized variants in increasing distance bins")
ggsave('Watershed_vs_distance.prioritized_or_not.png',p_compare,height = 8,width = 8)
```

```{r}
category_order <- c('coding variant','intron variant','1-1kb variant','1kb-10kb variant', '10-100kb variant')
p_compare <- revision_data_w_quantile %>% dplyr::select(nc_dist2gene, distance_categories, Watershed_quantile, CADDSV_quantile) %>% pivot_longer(cols = c(Watershed_quantile, CADDSV_quantile), 
               names_to = "Model", 
               values_to = "Score_quantile_bin") %>% mutate(Model = str_remove(Model, "_quantile$")) %>% 
  group_by(Model, distance_categories,Score_quantile_bin) %>% summarise(count = n()) %>%
  group_by(Model,distance_categories) %>% mutate(proportion = count / sum(count)) %>%
  ungroup() %>% 
  ggplot(aes(x = Score_quantile_bin, y = proportion, fill = factor(distance_categories,levels = category_order))) +
  geom_bar(stat = "identity", position = position_dodge(width = 0.9)) + 
  scale_fill_discrete() +
  facet_grid(rows = vars(Model)) +
  labs(x = "Score quantile bin", y = "fraction", fill = 'distance categories') +
  theme_minimal() +
  ggtitle("Grouped Bar Plot")
ggsave('Watershed_vs_CADD_distance.png',p_compare,height = 16,width = 16)
```

