---
title: "compare_starF_and_Arriba"
author: "bhaas"
date: '2024-07-20'
output: github_document
---

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


```{r}

starF_preds = read.csv("StarFusion_SGNEx_Illumina.tsv", header=T, sep="\t")

starF_preds %>% head()

```

```{r}

arriba_preds = read.csv("Arriba_SGNEx_Illumina.tsv", header=T, sep="\t")

arriba_preds %>% head()

```

```{r}

both_preds = full_join(starF_preds, arriba_preds, 
                       by=c('lex_ordered_fusion_name', 'sample', 'has_long_read_support', 'num_progs', 'progs'),
                       suffix=c('.starF', '.arriba')) %>%
    rename(LR_progs = progs, num_LR_progs = num_progs)

both_preds %>% head()

```


```{r}

# add indicator for what progs called it.
both_preds = both_preds %>% 
    mutate(progs = ifelse( (!is.na(`FusionName.starF`)) & (! is.na(`FusionName.arriba`)), 
                                                   "starF,arriba", "NA") ) %>%
    mutate(progs = ifelse(progs == "NA" & is.na(`FusionName.starF`), "arriba", progs) ) %>%
    
    mutate(progs = ifelse(progs == "NA" & is.na(`FusionName.arriba`), "starF", progs) )

both_preds

```

```{r}

both_preds %>% select(lex_ordered_fusion_name, progs, has_long_read_support) %>% unique() %>% group_by(progs, has_long_read_support) %>% tally()


```


```{r}
# compare arriba and starF regardless of LR support


both_preds %>% select(lex_ordered_fusion_name, progs) %>% unique() %>% 
    group_by(progs) %>% tally() %>% mutate(pct=prop.table(n)*100)

```



```{r}


# just those with long read support

both_preds %>% select(lex_ordered_fusion_name, progs, has_long_read_support) %>% unique() %>% 
    filter(has_long_read_support) %>%
    group_by(progs, has_long_read_support) %>% tally()  %>% ungroup() %>% mutate(pct=prop.table(n)*100)

```


```{r}
# how many with long read support in total?


both_preds %>% select(lex_ordered_fusion_name, progs, has_long_read_support) %>% unique() %>% 
    filter(has_long_read_support) %>% 
    select(lex_ordered_fusion_name) %>% unique() %>% nrow()


# 70
# but note, we recover another 3 with illumina support from our own K562 illumina sequencing., so 73 total.
```


```{r}

# just those WITHOUT long read support

both_preds %>% select(lex_ordered_fusion_name, progs, has_long_read_support) %>% unique() %>% 
    filter(! has_long_read_support) %>%
    group_by(progs, has_long_read_support) %>% tally() %>% ungroup() %>% mutate(pct=prop.table(n)*100)

```


```{r}

both_preds = both_preds %>% mutate(primary_fusion_name = ifelse(!is.na(`FusionName.starF`), `FusionName.starF`, `FusionName.arriba`))


```


```{r}

write.table(both_preds, file='SGNEx_Illumina_supported_fusions.tsv', quote=F, sep="\t", row.names=F)

```





