---
title: "UDN_nanopore_Watershed"
author: "Tanner Jensen"
date: "2023-11-26"
output: html_document
---

```{r setup, include=FALSE}
require(tidyverse)
require(magrittr)
library(arrow)
```
# Allele Frequency
```{r}
afs <- read_tsv("~/Downloads/SVAF_distribution_LongRead.tsv")
```

```{r}
afs
afs$rare <- afs$af < 0.0172
ggplot(afs, aes(af, fill=rare)) + geom_histogram() + theme_classic() + xlab("") + ylab("") + scale_fill_manual(values=c("grey", "#FFBC7C")) +
  theme(text=element_text(size=20), legend.position = "none")

```


```{r}
expression <- read_tsv("~/Downloads/EHMT1.tsv", col_names=c("sample", "gene", "zscore", "Batch", "UDN"))
expression$outlier <- abs(expression$zscore) > 2.5
ggplot(expression, aes(zscore, fill=outlier)) + geom_histogram() + theme_classic() + theme(legend.position = "none") + 
  scale_fill_manual(values=c("grey", "#2C70B0")) + ylim(0,150)
```

```{r colors}
SV_COLORS = c("DEL"="#C14E52","INS"="#2A2AA7","INV"="#EED54C","DUP"="#208179", "VNTR"="#AF2FD0")
SEQ_TECH = c("Short Read"="#A7CD59","Long Read"="#7F59CD")
OUTLIER_DIRECTION=list("over"="#FF8A75", "under"="#2EA3B8")
GENOME_ANNO_REGIOn=list("Noncoding"="#59C3AC","Coding"="#FF99AC")
```

```{r}
sv_full <- read_tsv("./LongVsShortStuff/SV_genotypes_by_rareness_and_length.shortread_updated.tsv")
sv_full_filt <- sv_full %>% filter(!(SVTYPE=="DUP"&length < 50)) %>% filter(!(SVTYPE=="INV"&length<50))
sv_counts <- sv_full %>% mutate(Category = ifelse(length < 50, yes="Large Indel",no="SV")) %>% group_by(Ind,Category,SeqTech,SVTYPE) %>% summarize(SV=n())
sv_counts_filt <- sv_full_filt %>% mutate(Category = ifelse(length < 50, yes="Large Indel",no="SV")) %>% group_by(Ind,Category,SeqTech,SVTYPE) %>% summarize(SV=n())
sv_counts %>% group_by(Ind, SeqTech, Category, SVTYPE) %>% summarize(SV = sum(SV)) %>% group_by(Category, SeqTech,SVTYPE) %>% summarise(median=median(SV))
sv_counts %>% group_by(Ind, SeqTech, SVTYPE) %>% summarize(SV=sum(SV)) %>% group_by(SVTYPE, SeqTech) %>% summarize(median=median(SV))
sv_counts_filt %>% group_by(Ind, SeqTech, SVTYPE) %>% summarize(SV=sum(SV)) %>% group_by(SVTYPE, SeqTech) %>% summarize(median=median(SV))
```

```{r}
sv_counts$SVTYPE %<>% factor(levels=c("DEL", "INS", "DUP", "INV"))
sv_counts$Category %<>% factor(levels=c("Large Indel", "SV"))
sv_counts$Category[sv_counts$SVTYPE=="DUP"] <- "SV"
sv_counts$Category[sv_counts$SVTYPE=="INV"] <- "SV"
sv_counts %>% group_by(Ind, SeqTech, Category, SVTYPE) %>% summarize(SV = sum(SV)) %>%
  ggplot(aes(SVTYPE, SV, fill=SeqTech)) + geom_boxplot() + scale_y_log10() + theme_classic() + scale_fill_manual(values=SEQ_TECH) + facet_grid(.~Category, scales="free", space="free_x") + ylab("SV count per individual") +
  theme(text=element_text(size=14), panel.grid.major.y=element_line()) + annotation_logticks(sides="l")
ggsave("./plots/figure1/SUP_FIG1.count_boxplot.pdf", width=10, height=7)

sv_counts_rare <- sv_full %>% filter(length > 50, rareness=="rare") %>% group_by(Ind, SeqTech,SVTYPE) %>% summarize(SV=n())
sv_counts_rare %>% group_by(Ind, SeqTech, SVTYPE) %>% summarize(SV = sum(SV)) %>% group_by(SVTYPE,SeqTech) %>% summarise(median=median(SV))
sv_counts %>% filter(rareness=="rare") %>% group_by(Ind, SeqTech, SVTYPE) %>% summarize(SV = sum(SV)) %>% group_by(SVTYPE,SeqTech) %>% summarise(median=median(SV))
sv_counts %>% filter(rareness=="rare") %>% group_by(Ind, SeqTech) %>% summarize(SV = sum(SV)) %>% group_by(SeqTech) %>% summarise(median=median(SV))
sv_counts %>% group_by(Ind, SeqTech, SVTYPE) %>% summarize(SV = sum(SV)) %>%
  ggplot(aes(SVTYPE, SV, fill=SeqTech)) + geom_boxplot() + scale_y_log10() + theme_classic() + scale_fill_manual(values=SEQ_TECH) + ylab("SV count per individual") +
  theme(text=element_text(size=14), panel.grid.major.y=element_line()) + annotation_logticks(sides="l")
ggsave("./plots/figure1/SUP_FIG1.count_boxplot.pdf", width=5, height=3)


rare_sv_counts <- fread("./LongVsShortStuff/SV_per_sample_stat_by_rareness.REVISION.tsv")
rare_sv_counts$SVTYPE %<>% factor(levels=c("DEL", "INS", "DUP", "INV"))
rare_sv_counts %>% filter(rareness=="rare") %>%
ggplot(aes(SVTYPE, SV, fill=SeqTech)) + geom_boxplot() + scale_y_log10() + theme_classic() + scale_fill_manual(values=SEQ_TECH) + 
  ylab("Rare variant count per individual") + xlab("Structural variant type") + theme(text=element_text(size=14), panel.grid.major.y = element_line()) + annotation_logticks(sides="l")
ggsave("./plots/figure1/fig1plotB.rare_boxplot.pdf", width=5, height=4)
```

```{r}
sv_lengths <- read_tsv("./LongVsShortStuff/SV_length_distribution_compare_long_and_short.short_read_updated.tsv")
sv_lengths %>% filter(length < 50) %>% group_by(cohort) %>% summarize(n())
sv_lengths %>% filter(length > 50) %>% group_by(cohort) %>% summarize(n())
77596/36985
120950/37517

table(sv_lengths$cohort)
sv_lengths %>% filter(SVTYPE=="INS") %>%
  ggplot(aes(length, color=cohort)) + geom_freqpoly(bins=60,position="identity") + 
  scale_x_log10(limits=c(30,NA), breaks=c(50,200,1e3,1e4,1e5), labels=c("50bp", "200bp", "1kb", "10kb", "100kb")) + 
  scale_y_log10(limits=c(1,NA)) + 
  theme_classic() + theme(text= element_text(size=14), panel.grid.major.y = element_line()) +
  scale_color_manual(values=SEQ_TECH) +
  geom_vline(xintercept = 50, linetype="dashed", color="black") +
  xlab("") + annotation_logticks(sides="l") + theme(legend.position = "none")
ggsave("./plots/figure1/fig1plotC.insertion.histogram.pdf", width=4, height=3)
```

```{r}
library("ggforce")
sv_lengths %>% filter(SVTYPE=="DEL") %>% 
  ggplot(aes(length, color=cohort)) + geom_freqpoly(bins=60,position="identity") + 
  theme_classic() + theme(text= element_text(size=14), legend.position = "none", panel.grid.major.y = element_line()) + 
  #scale_x_continuous(limits=c(25,1e8), breaks=c(50,200,1e3,1e4,1e5,1e6,1e7), labels=c("50bp","200bp","1kb","10kb", "100kb", "1Mb", "10Mb"), trans=trans_reverser('log10')) +
  scale_x_continuous(trans=trans_reverser('log10'),limits=c(1e7,30), breaks=c(5e6,1e6,1e5,1e4,1e3,200,50), labels=c("5Mb","1Mb","100Kb","10Kb","1Kb","200bp","50bp")) +
  scale_y_log10(position="right") + 
  geom_vline(xintercept = 50, linetype="dashed") +
  scale_color_manual(values=SEQ_TECH) + 
  xlab("") + annotation_logticks(sides="r")
ggsave('./plots/figure1/fig1plotC.deletion.histogram.pdf', width=5,height=3)
```

```{r}
sv_lengths %>% filter(SVTYPE=="DUP") %>%
  ggplot(aes(length, color=cohort)) + geom_freqpoly(bins=60,position="identity") + 
  scale_x_log10(limits=c(30,NA), breaks=c(100,1e3,1e4,1e5,1e6,1e7), labels=c("100bp", "1kb", "10kb", "100kb", "1Mb","10Mb")) + 
  scale_y_log10(limits=c(1,NA)) + 
  theme_classic() + theme(text= element_text(size=14), panel.grid.major.y = element_line()) +
  scale_color_manual(values=SEQ_TECH) + ggtitle("Duplications") +
  xlab("variant length") + annotation_logticks(sides="l")
ggsave("./plots/figure1/Dup.lengths.pdf", width=5,height=3)
sv_lengths %>% filter(SVTYPE=="INV") %>%
  ggplot(aes(length, color=cohort)) + geom_freqpoly(bins=60,position="identity") + 
  scale_x_log10(limits=c(30,NA), breaks=c(100,1e3,1e4,1e5,1e6,1e7), labels=c("100bp", "1kb", "10kb", "100kb", "1Mb","10Mb")) + 
  scale_y_log10(limits=c(1,NA)) + 
  theme_classic() + theme(text= element_text(size=14), panel.grid.major.y = element_line()) +
  scale_color_manual(values=SEQ_TECH) + ggtitle("Inversions") +
  xlab("variant length") + annotation_logticks(sides="l")
ggsave("./plots/figure1/Inv.lengths.pdf", width=5,height=3)
```

```{r}
sv_af <- read_tsv("./LongVsShortStuff/SVAF_distribution_LongRead.tsv")
min_af = round(min(sv_af$af),5)
sv_af %>%
  ggplot(aes(af, color=SVTYPE)) + geom_freqpoly(bins=60) + scale_y_log10() +
  scale_x_continuous(breaks=c(0.00075,.25,.5,.75,1.0), limits=c(0,1), labels=c(min_af, "0.25", "0.50", "0.75", "1.00")) +  
  theme_classic() + scale_color_manual(values=SV_COLORS) + xlab("SV Allele Frequency") +
  theme(panel.grid.major.y=element_line()) + annotation_logticks() 
ggsave('./plots/figure1/fig1plotB.freqpoly.pdf', width=4, height=3)
prop.table(table(sv_af$af < .01))
```

```{r}
short.sv_af <- read_tsv("./LongVsShortStuff/SVAF_distribution_ShortRead.tsv")
short.sv_af %>%
  ggplot(aes(af, color=SVTYPE)) + geom_freqpoly(bins=60) + scale_y_log10() + xlim(0,1) +
  theme_classic() + scale_color_manual(values=SV_COLORS) + xlab("SV Allele Frequency") +
  theme(panel.grid.major.y= element_line())
ggsave('./plots/figure1/SUP_FIG1.shortread.AF.freqpoly.pdf', width=5, height=4)
```

```{r}
sv_lengths %>% group_by(SVTYPE, cohort) %>% summarize(length=mean(length)) %>% filter(cohort=="Long Read")
```

```{r}
af_ascertainment <- read_tsv("./LongVsShortStuff/SV_af_by_svafotate_paragraph_OFP50.tsv")
af_ascertainment %>% group_by(SV) %>% summarize(SVAFotate=max(SVAFotate)) %>% pull(SVAFotate) %>% table %>% prop.table

head(af_ascertainment)
head(sv_full)
af_ascertainment
rare_var_counts_by_method <- rbind(
  af_ascertainment %>% filter(abs(SVLEN) > 50, SVAFotate_AF < 0.01, method!="Failed") %>% group_by(SUBJID,SVTYPE) %>% summarize(n = n()) %>% mutate(method="SVAFotate"),
  af_ascertainment %>% filter(abs(SVLEN) > 50, PARAGRAPH_AF < 0.01, method!="Failed") %>% group_by(SUBJID,SVTYPE) %>% summarize(n=n()) %>% mutate(method="PARAGRAPH 1000G"),
  sv_full %>% filter(abs(length) > 50, af < 0.01) %>% group_by(Ind, SVTYPE) %>% summarize(n=n()) %>% mutate(method="ADRC Jasmine-SV") %>% dplyr::rename("SUBJID"=Ind)
)
rare_var_counts_by_method %>% group_by(SUBJID, method) %>% summarize(n=sum(n)) %>% group_by(method) %>% summarize(n= median(n))
rare_var_counts_by_method$method %<>% factor(levels=c("SVAFotate", "PARAGRAPH 1000G", "ADRC Jasmine-SV"))
rare_var_counts_by_method$SVTYPE %<>% factor(levels=c("DEL","INS","DUP","INV"))
rare_var_counts_by_method %>%
  ggplot(aes(method, n, fill=SVTYPE)) +
  theme_classic() + theme(panel.grid.major.y=element_line(), panel.grid.minor.y=element_line(linewidth=.25)) +
  scale_fill_manual(values=SV_COLORS) +
  geom_boxplot() +  scale_y_log10()
ggsave("./plots/figure1/Supp_FIG1.Rare_variant_filtering.pdf", height=6, width=8)

```

## Enrichments
```{r}
enrich <- rbind(read_tsv("./UDNLongReadEnrichmentStuff/Result_matched_background_udn_enrichment_for_rareSV.blood.tsv"))
enrich
enrich %>% filter(`Z-threshold`==4) %>% mutate(CI=paste0(round(log_OR - 1.96*SE,3),"-",round(log_OR+1.96*SE,3)))
set.seed(0)
jit<-position_jitter(width=.06,height = 0)
ggplot(enrich, aes(`Z-threshold`, log_OR, ymin=log_OR-SE, ymax=log_OR+SE, color=SVTYPE,fill=SVTYPE)) + 
  geom_line(position = jit, linewidth=1) + 
  geom_pointrange(color="black",position=jit,size=.3,shape=21,fatten=20, alpha=.9) + 
  geom_hline(yintercept = 0, linetype="dashed") + 
  theme_classic() + 
  ylab("log(Odds Ratio)") +
  scale_color_manual(values=SV_COLORS) +
  scale_fill_manual(values=SV_COLORS) +
  guides(color = guide_legend(override.aes = list(size=1,linewidth=2))) +
  theme(legend.key.size = unit(.4,"inch"), panel.grid.major.y=element_line())
ggsave("./plots/Fig2.plotA.SVTYPE_enrichment.pdf", width=6,height=6)
```
```{r}
read_tsv("./UDNLongReadEnrichmentStuff/NEW_VNTR_enrichments.log_odds.txt") %>% filter(`Z-threshold` ==4) %>% mutate(CI= log_OR + 1.96*SE)
```

```{r}
overunder <- rbind(
  read_tsv("./UDNLongReadEnrichmentStuff/Result_matched_background_udn_underexpression_enrichment_for_rareSV.blood.tsv") %>% mutate(direction="under"),
  read_tsv("./UDNLongReadEnrichmentStuff/Result_matched_background_udn_overexpression_enrichment_for_rareSV.blood.tsv") %>% mutate(direction="over"),
  read_tsv("./UDNLongReadEnrichmentStuff/NEW_VNTR_enrichments.log_odds.txt") %>% filter(direction != "all")
)
overunder$SVTYPE %<>% factor(levels=c("DEL", "INS", "DUP","INV", "VNTR"))
overunder %>% filter(`Z-threshold`==4) %>%
  ggplot(aes(SVTYPE, log_OR, fill=direction,ymin=log_OR-SE,ymax=log_OR+SE)) +
    geom_pointrange(shape=21,fatten=10, position=position_dodge(width=.3)) +
  geom_hline(yintercept = 0, linetype="dashed") +
  scale_fill_manual(values=OUTLIER_DIRECTION) +
  theme_classic() + theme(panel.grid.major.y=element_line()) + 
  ylab("log(Odds Ratio)")
ggsave('./plots/Fig2.plotB.directional.SVTYPE_enrichment.pdf',width=4,height = 6)
overunder %>% filter(`Z-threshold`==4, direction=="over") %>% mutate(low=log_OR-1.96*SE, high=log_OR+1.96*SE)
```

```{r}
vntr_enrich <- read_tsv("./UDNLongReadEnrichmentStuff/NEW_VNTR_enrichments.log_odds.txt") %>% filter(direction != "all")
set.seed(0)
ggplot(vntr_enrich, aes(`Z-threshold`, log_OR, ymin=log_OR-SE, ymax=log_OR+SE, fill=direction)) + geom_line() + geom_pointrange(shape=21, position = position_jitter(width=.1), fatten=10) + 
  geom_hline(yintercept = 0, linetype="dashed")
ggsave("./plots/SUPP_FIG.VNTR.enrichment.pdf")

```


```{r}
short_vs_long <- read_tsv("LongVsShortStuff/Sample_matched_enrichment_result_Blood_long_vs_short.tsv")
short_vs_long <- read_tsv("LongVsShortStuff/Sample_matched_enrichment_result_Blood_long_vs_short.REVISION.tsv")
SEQ_TECH = c("short read"="#A7CD59","long read"="#7F59CD")
short_vs_long$SVTYPE %<>% factor(levels=c("DEL","INS","DUP","INV"))
short_vs_long %>%
  ggplot(aes(`Z-threshold`,log_OR,ymin=log_OR-SE,ymax=log_OR+SE,color=technology,fill=technology)) + geom_line() + 
  geom_pointrange(shape=21,fatten=5,color="black",position=position_jitter(width=.1)) +
  geom_hline(yintercept=0,linetype="dashed") + 
  ylim(-.7,NA) +
  facet_wrap(~SVTYPE,scales="free",ncol=1) + 
  theme_classic()  +
  scale_color_manual(values=SEQ_TECH) +
  scale_fill_manual(values=SEQ_TECH) 
ggsave('./plots/Fig2.plotC.short_vs_long.Enrichment.pdf',width=5,height=10)
```

```{r}
abc_element <- read_tsv("UDNLongReadEnrichmentStuff/UDN_ABC_enrichment_results.csv")
ABC_SVTYPE=c("has_ABC_rare_SV"="#59FFAC","has_coding_rare_SV"="#FF59AC")
abc_element %>% filter(tissue_enhancer=="All") %>%
  ggplot(aes(`Z-threshold`,log_OR,ymin=log_OR-SE,ymax=log_OR+SE,color=SVTYPE,fill=SVTYPE,linetype=tissue_expression)) +
  geom_line() +
  geom_pointrange(shape=21,color="black",linetype="solid", position=position_jitter(width=.08,height=0)) +
  geom_hline(yintercept = 0, linetype="dashed") +
  theme_classic()  +
  scale_fill_manual(values=ABC_SVTYPE)+
  scale_color_manual(values=ABC_SVTYPE)
ggsave('./plots/FIG2.plotD.ABC_enrichment.pdf')
#abc_element %>% filter(tissue_enhancer=="Fibroblast") %>%
  ggplot(aes(`Z-threshold`,log_OR,ymin=log_OR-SE,ymax=log_OR+SE,color=SVTYPE,fill=SVTYPE,linetype=tissue_expression)) +
  geom_line() +
  geom_pointrange(shape=21,color="black",linetype="solid", position=position_jitter(width=.08,height=0))
```

```{r}
tr_path_outliers <- read.csv("./TR_outlier_pathogenic_loci.csv")
tr_path_outliers
```

```{r}
ggplot(tr_path_outliers, aes(sample, RN, ymin=RN_lower, ymax=RN_upper, color=tool, shape=outlier)) + geom_pointrange(fatten=7,position=position_dodge2(width=1)) + facet_wrap(~Repeat, scale="free", ncol=4) + ylim(0,NA) + geom_hline(aes(yintercept=pathogenic_threshold), linetype="dashed", color="red") + scale_color_manual(values=c("#A7CD59","#7F59CD")) + theme_classic() + 
  theme(panel.grid.major.y = element_line(),
        axis.text.x = element_text(hjust=1,vjust=1,angle=45)) + xlab("") + ylab("repeat copy number")
ggsave('./TR_path_outliers.pdf', height=10, width=8)
```

```{r}
library(ggupset)
library(cowplot)
library(dplyr)
library(tidyr)
library(tibble)
library(stringr)
library(ggplot2)
library(ggarange)
```



```{r}
prioritization_df <- read_tsv("./Fig4Stuff/Lolipop_raw_data.tab.tsv")
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)))
e
# Stripcharts with summary statistics
# Change color by dose groups
custom_colors <- c("WGS-only" = "#FFBC7B", "WGS-RNA" = "#1E88E5", "WGS-HPO" = "#004D40", "WGS-RNA-HPO" = "#D81B60")
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() + 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') +
  theme_classic() + 
  theme(panel.grid.major.x=element_line(), panel.grid.minor.x=element_line(linewidth=.2), legend.position="bottom") +
  guides(colour=guide_legend(nrow=1))
e1

e2 <- ggplot(lolipop_plot_data, aes(y = factor(filter,levels = orders))) + geom_bar(fill="grey") + labs(x = 'Patient Count') + theme_classic() + theme(aspect.ratio=3/1,
                                                                                                                          panel.grid.major.x=element_line(),
                                                                                                                          axis.text.y=element_blank(),
                                                                                                                          axis.title.y=element_blank(),
                                                                                                                          axis.ticks.y=element_blank())
e2
e1
ggsave("./plots/Fig4A.SV-gene_prioritization.pdf", height=10, width=10)
e2
ggsave("./plots/Fig4A.prioritized_patient_count.pdf", width=3, height=9) 
e2
figure5a <- ggarrange(e1,e2,nrow = 1,widths = c(10,3),heights = 10)
figure5a
figure5a %>% ggsave(file="Figure4/filter_plot.svg", width=11, height=5, dpi=300)

```

```{r}
lr_sr_df <- read_delim('./Fig4Stuff/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('./Fig4Stuff/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()
prioritization_df_compare_sr_lr <- prioritization_df %>% mutate(unique_to_long_read = ifelse(SV %in% lr_sr_shared_updated, F, T))
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(fill="grey") +
    theme_classic() +
  theme(panel.grid.major.y=element_line(), panel.grid.minor.y=element_line(linewidth=.25)) +
    scale_x_upset()  + scale_y_continuous(position = "right")
upset_long_vs_short

upset_long_vs_short %>% ggsave(file="./plots/Fig4B.Upset_watershed_caddsv_lr_sr.pdf", width=5, height=4)
upset_long_vs_short
```

```{r}
proportion_plot <- read_tsv("./LongVsShortStuff/SR_LR_proportion_plot_df.40.REVISION.tsv")
colnames(proportion_plot) <- c("SVTYPE", "LR_vs_SR", "Percentage")
proportion_plot$LR_vs_SR %<>% factor(levels=c("LR_only", "both", "SR_only"))
proportion_plot$SVTYPE %<>% factor(levels=c("DEL", "INS", "DUP", "INV"))
proportion_plot %>% 
  ggplot(aes(SVTYPE, Percentage, label=paste0(round(Percentage*100,1),"%"), fill=LR_vs_SR)) + geom_col() + scale_fill_manual(values=c(LR_only="#7F59CD", SR_only="#A7CD59", both="#F58220")) +
  geom_text(size=5, position = position_stack(vjust=0.5), color="white") 
ggsave("./plots/Supp_FIG.SR_vs_LR_outlier_variants.pdf",height=9, width=8)
```


```{r}
svs <- read_tsv("~/Downloads/ALNG.GRCh38.consensus_SVs.jasmine_irisRefined.AF_annotated.txt")
svs.formatted <- svs %>% select(seqnames, start, end, type, size, ALNG_AF) %>% dplyr::rename(chrom="seqnames", length="size", AF="ALNG_AF")
write.table(svs.formatted, "~/Documents/Research/UDN_Watershed_Figures/Supp_Table4_SV_coordinates_AFs.txt", sep="\t", col.names=T, row.names=F, quote=F)
```

```{r}
trs <- read_tsv("~/Downloads/combined_TR_mean_variance_Technology.table.txt")
trs %>% filter(Technology!= "Long Reads (vamos)")
trs %<>% mutate(ID = ifelse(is.na(VNTR_ID), name, VNTR_ID))
colnames(trs)

nrow(trs)
trs_formatted <- trs %>% select(ID, RS, Technology, bin,median_RL, mean_RL, var_RL)
write.table(trs_formatted, file="~/Documents/Research/UDN_Watershed_Figures/Supp_Table3.Tandem_repeat_copy_number_variance_mean.txt", sep="\t", col.names=T, row.names=F, quote=F)
table(trs_formatted$Technology)
trs_formatted
```


