---
title: "ExamineFusionPredCounts"
author: "bhaas"
date: '2024-07-21'
output: github_document
---

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


```{r}

PROGS = c('ctat-LR-fusion', 'JAFFAL', 'LongGF', 'fusionseeker', 'pbfusion')

```

```{r}

fusion_preds = read.csv("__bmark_min-1-read/data/preds.collected.gencode_mapped.wAnnot.gz", header=T, sep="\t") %>% 
    select(sample, prog, num_reads, fusion)

fusion_preds$prog = factor(fusion_preds$prog, levels=PROGS)

fusion_preds =  fusion_preds %>% rowwise() %>% mutate(lex_ordered_fusion_name = paste0(sample, "|", paste0(collapse="--", sort(str_split(fusion, "--")[[1]])))) 

```


```{r}

fusion_preds = fusion_preds %>% rowwise() %>% mutate(lex_fusion = paste0(collapse="--", sort(str_split(fusion, "--")[[1]])))

```


```{r}

# remove messy fusions
# defined as having at least 3 reads and predicted across mult samples

messy_fusions = fusion_preds %>% filter(num_reads >= 3) %>% select(sample, lex_fusion) %>% unique() %>%
    group_by(lex_fusion) %>% tally() %>% filter(n>1) %>% pull(lex_fusion)

fusion_preds = fusion_preds %>% filter(! lex_fusion %in% messy_fusions)    

length(messy_fusions)

```





```{r}
# arriba fusion preds
arriba_fusion_preds = read.csv("3b.1.IlluminaTruSeqDepMap9Lines/all_arriba_high_preds.tsv", sep="\t", header=T, stringsAsFactors = F)
arriba_fusion_preds %>% select(lex_ordered_fusion_name) %>% unique() %>% nrow()

# 282 arriba fusions
```


```{r}
# starF preds
starF_preds = read.csv("3b.1.IlluminaTruSeqDepMap9Lines/all_STARF_preds.tsv", sep="\t", header=T, stringsAsFactors = F)
starF_preds %>% select(lex_ordered_fusion_name) %>% unique() %>% nrow()


# 266 StarF preds

```



```{r}

illumina_preds = bind_rows(arriba_fusion_preds %>% select(sample, lex_ordered_fusion_name) %>% unique() %>% mutate(prog='arriba'),
                           starF_preds %>% select(sample, lex_ordered_fusion_name) %>% unique() %>% mutate(prog='starF')
                           )

illumina_preds = illumina_preds %>% group_by(sample, lex_ordered_fusion_name) %>% mutate(progs = paste0(collapse=',', prog)) %>%
    ungroup() %>% select(sample, lex_ordered_fusion_name, progs) %>% unique()

illumina_preds %>% head()

```

```{r}

table(illumina_preds$progs)

# 190 +  92 +  174  = 456

```


```{r}

illumina_preds %>% group_by(sample, progs) %>% tally() 

```




```{r}

illumina_preds %>% group_by(sample, progs) %>% tally() %>%
    ggplot(aes(x=sample, y=n, fill=progs)) + geom_col()


```


```{r}


fusion_preds = fusion_preds %>% mutate(illumina_supported = lex_ordered_fusion_name %in% illumina_preds$lex_ordered_fusion_name)
nrow(fusion_preds)
# take top supported fusion call per prog in case mult isoforms
fusion_preds = fusion_preds %>% group_by(prog, lex_ordered_fusion_name) %>% arrange(desc(num_reads)) %>% filter(row_number() == 1) %>% ungroup()

nrow(fusion_preds)
```


```{r}

fusion_counts_df = NULL

for (i in seq(1,10)) {
    
    fusion_counts = fusion_preds %>% filter(num_reads >= i) %>% group_by(prog, illumina_supported) %>% tally() %>% mutate(min_reads=i)

    fusion_counts_df = bind_rows(fusion_counts_df, fusion_counts)
    
}

fusion_counts_df %>% head()

```

```{r}

TP_vs_FP_illum_supported_barplot = fusion_counts_df %>% ggplot(aes(x=prog, y=n)) + geom_col(aes(fill=illumina_supported)) +
    theme_bw() +
    facet_wrap(~min_reads, scale='free_y') +
    theme(axis.text.x = element_text(angle = 90, hjust = 1)) 

TP_vs_FP_illum_supported_barplot
```
```{r}

ggsave(TP_vs_FP_illum_supported_barplot, filename="TP_vs_FP_illum_supported_barplot.svg", width=7, height=4.5)

```

```{r}

frac_fusion_illum_supported = fusion_counts_df %>% spread(key=illumina_supported, value=n) %>% 
    mutate(frac_illum_supported = `TRUE` / (`TRUE`+`FALSE`)) 


```


```{r}

frac_supported_vs_min_reads_plot = frac_fusion_illum_supported %>% ggplot(aes(x=min_reads, y=frac_illum_supported)) + geom_line(aes(groups=prog, color=prog)) +
    geom_point(aes(color=prog)) +
    theme_bw()

frac_supported_vs_min_reads_plot
```



```{r}

#ggsave(frac_supported_vs_min_reads_plot, filename="frac_supported_vs_min_reads_plot.svg", width=7, height=4.5)


```


# examine min 2 agree profile

```{r}

fusion_preds_min_agree = fusion_preds %>% select(prog, lex_ordered_fusion_name) %>% unique() %>%
    group_by(lex_ordered_fusion_name) %>% tally() %>% filter(n>1)

fusion_preds_min_agree %>% head()

```

```{r}

fusion_preds_min2_agree_names = fusion_preds %>% mutate(min2_agree = lex_ordered_fusion_name %in% (fusion_preds_min_agree %>% filter(n>=2) %>% pull(lex_ordered_fusion_name)) )


fusion_preds_min3_agree_names = fusion_preds %>% mutate(min3_agree = lex_ordered_fusion_name %in% (fusion_preds_min_agree %>% filter(n>=3) %>% pull(lex_ordered_fusion_name)) )



```

```{r}

fusion_counts_min_agree_df = NULL

for (i in seq(1,10)) {
    
    for (min_progs_agree in seq(2,4)) {
        
         fusions_with_min_reads = fusion_preds %>% filter(num_reads >= i) 
            num_fusions_min_reads = fusions_with_min_reads %>% select(lex_ordered_fusion_name) %>% unique() %>% nrow()
    
        fusion_counts_min_agree = fusion_preds %>% filter(num_reads >= i) %>% 
             group_by(lex_ordered_fusion_name) %>% tally() %>% filter(n>=min_progs_agree) %>% nrow() 

            fusion_counts_min_agree_df = bind_rows(fusion_counts_min_agree_df, 
                                            data.frame(min_read_threshold = i,
                                                       min_progs_agree = min_progs_agree,
                                                       num_fusions_min_reads = num_fusions_min_reads,
                                                       num_fusions_min_agree = fusion_counts_min_agree,
                                                       frac_fusions_min_agree = fusion_counts_min_agree/num_fusions_min_reads
                                                        ) 
                                            
         )
   }   
}

fusion_counts_min_agree_df 

```


```{r}

# fraction of all predictions that are agreed upon given minimum read thresholds applied

frac_fusions_min_progs_agree_plot = fusion_counts_min_agree_df  %>%
    mutate(min_progs_agree = factor(min_progs_agree)) %>%
    ggplot(aes(x=factor(min_read_threshold), y=frac_fusions_min_agree)) + 
    theme_bw() +
    geom_point(aes(color=min_progs_agree)) + geom_line(aes(group=min_progs_agree, color=min_progs_agree))

frac_fusions_min_progs_agree_plot

```

key observations:
- increasing read support yields increasingly agreeing predictions across programs with the exception of min 1 read support for min-2 and min-3 progs agreeing.
- increasing hte number of programs that must agree reduces the fraction that agree at the min-1 read evidence threshold point.





```{r}

ggsave(frac_fusions_min_progs_agree_plot, filename="frac_fusions_min_progs_agree_plot.svg", width=7, height=4.5)


```

# examine truth sets overall fraction


```{r}

fusion_counts_min_agree_vs_illum_supported = NULL

for (i in seq(1,10)) {
    
    for (min_progs_agree in seq(2,4)) {
        
        num_fusion_preds = fusion_preds %>% filter(num_reads >= i) %>%
            group_by(lex_ordered_fusion_name, illumina_supported) %>% tally() %>% filter(n > min_progs_agree) %>% ungroup() %>%
            group_by(illumina_supported) %>% tally(name='num_fusions')
        
        fusion_counts_min_agree_vs_illum_supported = bind_rows(fusion_counts_min_agree_vs_illum_supported, 
                                                           num_fusion_preds %>% mutate(min_read_support = i, min_progs_agree = min_progs_agree)
                                                           )
    }
}

fusion_counts_min_agree_vs_illum_supported

```



```{r}

frac_illum_supported_by_min_progs_agree_plot = fusion_counts_min_agree_vs_illum_supported %>% 
    spread(key=illumina_supported, value=num_fusions) %>% mutate(frac_illum_supported = `TRUE`/(`TRUE`+`FALSE`)) %>%
    ggplot(aes(x=factor(min_read_support), y=frac_illum_supported)) + theme_bw() +
    geom_point(aes(color=factor(min_progs_agree))) + 
    geom_line(aes(group=min_progs_agree, color=factor(min_progs_agree)))
   

frac_illum_supported_by_min_progs_agree_plot

```


Increasing the minimum read support before defining proxy truth sets increases the fraction of fusions supported by orthogonal illumina read evidence.

Increasing the number of programs that must agree, particularly at 3 or more reads as evidence, increases the relative illumina-support among the proxy truth set.


```{r}

ggsave(frac_illum_supported_by_min_progs_agree_plot, filename="frac_illum_supported_by_min_progs_agree_plot.svg", width=7, height=4.5)

```




