---
title: "SGNex_ONT_eval"
author: "bhaas"
date: '2024-07-24'
output: github_document
---

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


```{r}

PREDS_FILE_ALL = "data/preds.collected.gencode_mapped.wAnnot.filt.pass.proxy_assignments.gz"

# or use the version with ctat-LRF restricted to spliced fusions for comparison.
PREDS_FILE_fuzzy_restricted = "data/preds.collected.gencode_mapped.wAnnot.filt.pass.proxy_assignments.FUZZY_RESTRICTED.gz"

```

Proxy fusion names

```{r}

parse_preds_file = function(PREDS_FILE) {

    read_support = read.csv(PREDS_FILE, header=T, sep="\t")
    
    read_support$sample = str_replace(read_support$sample, "-EV_directRNA", "_EVdirectRNA")
    
    return(read_support)    
}


```

```{r}

regular_run_data = parse_preds_file(PREDS_FILE_ALL)

regular_run_data %>% head()


```

```{r}

fuzzy_restricted_run_data = parse_preds_file(PREDS_FILE_fuzzy_restricted)


fuzzy_restricted_run_data %>% head()

```

```{r}

regular_run_data$runtype = "regular"
fuzzy_restricted_run_data$runtype  = "fuzzy_restricted"

read_support = bind_rows(regular_run_data %>% mutate(runtype = "regular"),
                         fuzzy_restricted_run_data %>% mutate(runtype = "fuzzy_brkpt_restricted") )


```

```{r}

read_support = read_support %>% select(-c(mapped_gencode_A_gene_list, mapped_gencode_B_gene_list))

```



```{r}

# use lex-sorted fusion name

make_lex_sorted_fusion_name = function(fusion_name) {
    
    sample_name = str_split(fusion_name, "\\|")[[1]][1]
    fusion_pair = str_split(fusion_name, "\\|")[[1]][2]
    
    sorted_gene_pair = sort(str_split(fusion_pair,"--")[[1]])
    
    lex_sorted_fusion_name = paste0(sample_name, "|", sorted_gene_pair[1], "--", sorted_gene_pair[2])
    
    return(lex_sorted_fusion_name)
}

read_support = read_support %>% rowwise() %>% mutate(lex_sorted_fusion_name = make_lex_sorted_fusion_name(proxy_fusion_name))

```


```{r}

# restrict to dominant isoform in each sample

nrow(read_support)

read_support_sample = read_support %>% group_by(lex_sorted_fusion_name, prog, sample, runtype) %>% arrange(desc(num_reads)) %>%
    filter(row_number() == 1) %>%
    ungroup()

nrow(read_support_sample)
```



```{r}

validated_fusions = data.frame(validated_fusion = c(
"MCF7|BCAS4--BCAS3",
"MCF7|ARFGEF2--SULF2",
"MCF7|RPS6KB1--VMP1",
"MCF7|GCN1L1--MSI1",
"MCF7|AC099850.1--VMP1",
"MCF7|SMARCA4--CARM1",
"MCF7|SLC25A24--NBPF6",
"MCF7|USP31--CRYL1",
"MCF7|TBL1XR1--RGS17",
"MCF7|TAF4--BRIP1",
"MCF7|RPS6KB1--DIAPH3",
"MCF7|AHCYL1--RAD51C",
"MCF7|ABCA5--PPP4R1L",
"MCF7|C16orf45--ABCC1",
"MCF7|C16orf62--IQCK",
"MCF7|TXLNG--SYAP1",
"MCF7|MYO6--SENP6",
"MCF7|POP1--MATN2",
"MCF7|GATAD2B--NUP210L",
"MCF7|ESR1--CCDC170",
"MCF7|DEPDC1B--ELOVL7",
"MCF7|ATXN7L3--FAM171A2",
"MCF7|SYTL2--PICALM",
"MCF7|ADAMTS19--SLC27A6",
"MCF7|ARHGAP19--DRG1",
"MCF7|MYO9B--FCHO1",
"MCF7|PAPOLA--AK7",
"MCF7|ANKS1A--UHRF1BP1",
"MCF7|ATP1A1--ZFP64",
"MCF7|B3GNTL1--SLC9A8",
"MCF7|BCAS3--AMPD1",
"MCF7|BCAS3--ATXN7",
"MCF7|BCAS4--ZMYND8",
"MCF7|CHEK2--XBP1",
"MCF7|KCND3--PPM1E",
"MCF7|RP11-145E5.5--CDKN2B-AS1",
"MCF7|NAV1--GPR37L1",
"MCF7|NCOA3--SULF2",
"MCF7|PLCG1--TOP1",
"MCF7|PNPLA7--DPH7",
"MCF7|PTPRG--CCDC129",
"MCF7|PTPRN2--FAM62B",
"MCF7|RAD51C--ATXN7",
"MCF7|RSBN1--AP4B1-AS1",
"MCF7|SGPP2--ULK4",
"MCF7|SULF2--PRICKLE2",
"MCF7|TEX14--PTPRG",
"MCF7|TOP1--CR593014",
"MCF7|TSPAN9--TEAD4",
"MCF7|UBE2V1--TBX2",
"MCF7|VAV3--AP4B1-AS1",
"MCF7|ZMYND8--USP32",
"MCF7|MYH9--EIF3D",

"K562|NUP214--XKR3",
"K562|PRIM1--NACA",
"K562|BRK1--VHL",
"K562|ACCS--EXT2",
"K562|SLC29A1--HSP90AB1",
"K562|BCR--ABL1",

"A549|SCAMP2--WDR72",
"A549|RICTOR--FYB"

) )


validated_fusions = validated_fusions %>% rowwise() %>% mutate(lex_sorted_validated_fusion = make_lex_sorted_fusion_name(validated_fusion))

```


```{r}

lex_sorted_validated_fusions = validated_fusions %>% select(lex_sorted_validated_fusion) %>% unique() %>% pull(lex_sorted_validated_fusion)


```


```{r}

# incorporate the validated info
read_support_sample = read_support_sample %>% rowwise() %>% 
    mutate(validated_fusion = (lex_sorted_fusion_name %in% lex_sorted_validated_fusions) )


```


```{r}

read_support_sample %>% filter(validated_fusion) %>% 
    select(sample, lex_sorted_fusion_name, prog) %>% unique() %>%
    group_by(lex_sorted_fusion_name) %>% 
    arrange(prog) %>% 
    mutate(progs = paste(prog, collapse=","), num_progs = n()) %>% 
    select(sample, lex_sorted_fusion_name, progs, num_progs) %>%
    unique() %>%
    arrange(sample, desc(num_progs))

```



# incorporate illumina support info

```{r}

illumina_support_info = read.csv("6.1.SGNex_Matched_Illumina/SGNEx_Illumina_supported_fusions.tsv", header=T, sep="\t")

illumina_support_info %>% head()

```

```{r}

read_support_sample = left_join(read_support_sample,
                             illumina_support_info %>% select(lex_ordered_fusion_name, progs) %>% 
                                 unique() %>%
                                 rename(matched_illumina = progs),
                             by=c('lex_sorted_fusion_name'='lex_ordered_fusion_name') )

```

```{r}

read_support_sample %>% filter(validated_fusion | !is.na(matched_illumina)) %>% 
    select(sample, lex_sorted_fusion_name, validated_fusion, matched_illumina, prog) %>% unique() %>%
    group_by(lex_sorted_fusion_name) %>% 
    arrange(prog) %>% 
    mutate(progs = paste(prog, collapse=","), num_progs = n()) %>% 
    select(sample, lex_sorted_fusion_name, validated_fusion, matched_illumina, progs, num_progs) %>%
    unique() %>%
    arrange(sample, desc(num_progs))

```


# incorporate our separate Illumina supported fusions

```{r}

our_illumina_reads = read.csv("../3.DepMap9Lines/3b.DepMap9Lines_Benchmarking/3b.1.IlluminaTruSeqDepMap9Lines/Illumina_supported_fusions.tsv", header=T, sep="\t")


read_support_sample = left_join(read_support_sample, 
                             our_illumina_reads %>% select(lex_ordered_fusion_name, progs) %>% 
                                 unique() %>%
                                 rename(other_illumina = progs),
                             by=c('lex_sorted_fusion_name'='lex_ordered_fusion_name') )
```


```{r}

read_support_sample %>% filter(validated_fusion | !is.na(matched_illumina) | !is.na(other_illumina) ) %>% 
    select(sample, lex_sorted_fusion_name, validated_fusion, matched_illumina, other_illumina, prog) %>% unique() %>%
    group_by(lex_sorted_fusion_name) %>% 
    arrange(prog) %>% 
    mutate(progs = paste(prog, collapse=","), num_progs = n()) %>% 
    select(sample, lex_sorted_fusion_name, validated_fusion, matched_illumina, other_illumina, progs, num_progs) %>%
    unique() %>%
    arrange(sample, desc(num_progs))

```





```{r}

# which validated fusions do we find?

read_support_sample %>% filter(validated_fusion) %>% select(lex_sorted_fusion_name) %>% unique()

```




```{r}

# which validated fusions are missing?


setdiff(lex_sorted_validated_fusions, read_support_sample %>% select(lex_sorted_fusion_name) %>% unique() %>% pull(lex_sorted_fusion_name))


```


```{r}

read_support_sample = read_support_sample %>% mutate(as_truth = (validated_fusion | !is.na(matched_illumina) | !is.na(other_illumina)) )

```


```{r}
read_support_sample %>% filter(as_truth) %>% 
    select(lex_sorted_fusion_name, validated_fusion, matched_illumina, other_illumina) %>% 
    unique()


# 79 'truth' fusions
```

# score fusions

```{r}

# score fusions

read_support_sample_scored = read_support_sample %>% 
     group_by(prog, runtype) %>%
     arrange(desc(num_reads), desc(as_truth)) %>% 
    mutate(num_valid = cumsum(as_truth), num_other = cumsum(! as_truth)) %>%
    ungroup()

```

```{r}
read_support_sample_scored %>% filter(prog == 'ctat-LR-fusion')  %>% head(50)
```


```{r}

read_support_sample_scored$prog = factor(read_support_sample_scored$prog, levels=c('ctat-LR-fusion',
                                                                             'JAFFAL',
                                                                             'LongGF',
                                                                             'fusionseeker',
                                                                             'pbfusion'
                                                                             ))
```


```{r}

alpha_val = 0.2

SGNEx_fusion_plot = read_support_sample_scored %>% 
    select(prog, num_other, num_valid, runtype) %>% 
    ggplot(aes(x=num_other, y=num_valid)) + 
    geom_point(aes(color=prog), alpha=alpha_val) + 
    geom_line(aes(color=prog)) +
    scale_x_continuous(trans='log10') + theme_bw() +
    facet_wrap(~runtype)

SGNEx_fusion_plot
```


```{r}
# simplify plot for upload into inkscape

num_valid_vs_num_other = read_support_sample_scored %>% 
    select(prog, num_other, num_valid, runtype) 


num_valid_vs_num_other_lte_10k = num_valid_vs_num_other %>% group_by(runtype) %>% filter(num_other <= 10000) %>% ungroup()
num_valid_vs_num_other_gt_10k = num_valid_vs_num_other %>% group_by(runtype) %>% filter(num_other > 10000) %>% ungroup()



```

```{r}


num_valid_vs_num_other_gt_10k  = num_valid_vs_num_other_gt_10k  %>% group_by(runtype) %>% sample_n(1000) %>% ungroup()

```

```{r}
num_valid_vs_num_other  = bind_rows(num_valid_vs_num_other_lte_10k, num_valid_vs_num_other_gt_10k)

SGNEx_fusion_plot = num_valid_vs_num_other %>% 
    select(prog, num_other, num_valid, runtype)  %>% 
    
    mutate(runtype = factor(runtype, levels=c('regular', 'fuzzy_brkpt_restricted'))) %>%
    
    ggplot(aes(x=num_other, y=num_valid)) + 
    geom_point(aes(color=prog), alpha=alpha_val) + 
    geom_line(aes(color=prog)) +
    scale_x_continuous(trans='log10') + theme_bw() + facet_wrap(~runtype)

SGNEx_fusion_plot

```




```{r}

ggsave(SGNEx_fusion_plot, filename="SGNEx_fusion_plot.svg", width=7, height=3.5)


```




```{r}

# examine Venn for valid fusions

read_support_sample_scored  %>% filter(as_truth) %>% group_by(lex_sorted_fusion_name, runtype) %>%
    arrange(prog) %>%
    mutate(progs = paste(collapse=",", prog)) %>%
    ungroup() %>%
    select(runtype, lex_sorted_fusion_name, progs) %>% unique() %>% 
    arrange(runtype, lex_sorted_fusion_name, progs)


```

```{r}

read_support_sample_scored  %>% filter(as_truth) %>% group_by(lex_sorted_fusion_name, runtype) %>%
    arrange(prog) %>%
    mutate(progs = paste(collapse=",", prog)) %>%
    ungroup() %>%
    select(lex_sorted_fusion_name, progs, runtype) %>% unique() %>% 
    group_by(progs, runtype) %>% tally() %>%
    arrange(desc(n))
    


```







```{r}

read_support_sample_scored  %>% filter(as_truth) %>% group_by(lex_sorted_fusion_name) %>%
    arrange(prog, runtype) %>%
    mutate(progs = paste(collapse=",", prog)) %>%
    ungroup() %>%
    select(lex_sorted_fusion_name, progs, runtype) %>% unique() %>% 
    group_by(progs, runtype) %>% tally() %>%
    ggplot(aes(x=progs, y=n)) + 
    geom_col() +
     theme(axis.text.x = element_text(angle = 90, hjust = 1)) +
    facet_wrap(~runtype) +
    ylab("# as-truth fusions")

```


```{r}
# get counts of valid fusions

read_support_sample_scored  %>% filter(validated_fusion) %>% 
    select(lex_sorted_fusion_name, prog, runtype) %>% unique() %>% group_by(prog, runtype) %>% tally() %>%
    arrange(desc(n))
    
```

```{r}

read_support_sample %>% ungroup() %>% select(sample) %>% unique()


```

```{r}

# total number of fusion preds

read_support_sample %>% ungroup() %>% select(prog, lex_sorted_fusion_name, runtype) %>% 
    unique() %>% group_by(prog, runtype) %>% tally() %>% 
    arrange(desc(n))


```

```{r}

# write supp data table

write.table(read_support_sample_scored %>% 
                filter(runtype == "regular") %>%
                arrange(prog, desc(num_reads), num_valid, num_other) %>%
                select(proxy_fusion_name, sample, prog, num_reads, validated_fusion, matched_illumina, other_illumina), 
            file="Table_Sx-ONT_fusions_by_method.tsv", 
            sep="\t", quote=F, row.names=F)

```



```{r}

lex_sorted_validated_fusions = read_support_sample_scored  %>% filter(as_truth) %>% 
    select(lex_sorted_fusion_name, validated_fusion, matched_illumina, other_illumina ) %>% 
    unique() %>%
    arrange(lex_sorted_fusion_name)

write.table(lex_sorted_validated_fusions, 
            file='SGNEx-as_truth_fusions.lex_ordered.tsv', quote=F, row.names=F, col.names=T, sep="\t")

lex_sorted_validated_fusions

```


```{r}

lex_sorted_validated_fusions %>%
    mutate(has_matched_illumina = ! is.na(matched_illumina),
           has_other_illumina = ! is.na(other_illumina)) %>%
    group_by(validated_fusion, has_matched_illumina, has_other_illumina) %>% tally()


```

```{r}
# how many validated fusions included?

lex_sorted_validated_fusions %>% filter(validated_fusion) %>% nrow()

# 29 validated fusions

# 75-29 = 46 supplemented fusions w/ illumina support

# 3 of the illumina supported fusions derive from our own K562 illumina

```


```{r}

# how many with matched illumina

lex_sorted_validated_fusions %>% filter(! is.na(matched_illumina)) %>% nrow()


```


```{r}

lex_sorted_validated_fusions %>%  filter(! is.na(matched_illumina)) %>% group_by(matched_illumina) %>% tally()


```


```{r}

# 73 fusiosn w/ Illumina support

lex_sorted_validated_fusions %>%  filter(! (is.na(matched_illumina) & is.na(other_illumina))) %>% nrow()

```


```{r}

lex_sorted_validated_fusions %>%  filter(! (is.na(matched_illumina) & is.na(other_illumina))) %>% group_by(matched_illumina, other_illumina) %>% tally()

```

# counting trusted fusions by method

```{r}

read_support_sample %>% group_by(prog, runtype) %>% filter(as_truth) %>% tally() %>% arrange(runtype, desc(n))


```

# count the validated ones

```{r}

read_support_sample %>% group_by(prog, runtype) %>% filter(validated_fusion) %>% tally() %>% arrange(runtype, desc(n))


```


```{r}

trusted_fusions_found = read_support_sample %>%
    filter(as_truth) %>%
    select(prog, proxy_fusion_name, validated_fusion, matched_illumina, other_illumina)  %>%
    unique() %>%
    group_by(proxy_fusion_name, validated_fusion, matched_illumina, other_illumina) %>%
    arrange(prog) %>%
    mutate(progs = paste(collapse=",", prog)) %>%
    select(proxy_fusion_name, validated_fusion, matched_illumina, other_illumina, progs) %>% unique() 
    
trusted_fusions_found

```

```{r}

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


```




```{r}

read_support_sample_scored %>% filter(prog == "ctat-LR-fusion") %>%
    filter(runtype == "regular") %>%
    filter(! as_truth) %>% 
    select(proxy_fusion_name, fusion, num_reads)



```




```{r}
# Examine Venn for valid and illumina support for trusted fusions

trusted_fusions_found %>% mutate(matched_illumina = ifelse(is.na(matched_illumina), other_illumina, matched_illumina)) %>%
    group_by(validated_fusion, matched_illumina) %>% tally()

# Venn for Figure 6a

```

```{r}

# include lex sorted fusion name
trusted_fusions_found_adj = trusted_fusions_found %>% mutate(matched_illumina = ifelse(is.na(matched_illumina), other_illumina, matched_illumina)) %>%
    rowwise() %>% mutate(lex_ordered_fusion_name = make_lex_sorted_fusion_name(proxy_fusion_name))


trusted_fusions_found_adj  = left_join(trusted_fusions_found_adj,
                                       illumina_support_info %>% select(lex_ordered_fusion_name, FFPM) %>% 
                                           group_by(lex_ordered_fusion_name) %>% arrange(desc(FFPM)) %>% filter(row_number() == 1) %>% ungroup(),
                                       by='lex_ordered_fusion_name')


```


```{r}

trusted_fusions_found_adj %>% 
    filter(grepl("starF", matched_illumina)) %>%
    ggplot(aes(x=matched_illumina, y=FFPM)) + geom_jitter(aes(color=matched_illumina)) +
    theme_bw() +
    scale_y_continuous(trans='log10')


```



```{r}


wilcox.test(trusted_fusions_found_adj %>% filter(matched_illumina == "starF,arriba") %>% pull(FFPM),
            trusted_fusions_found_adj %>% filter(matched_illumina == "starF") %>% pull(FFPM),
            alternative = 'greater')


```







