---
title: "CTAT_DepMap9Lines"
author: "bhaas"
date: '2024-02-01'
output: github_document
---

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


```{r}

files = paste0("data/", list.files("data/", "*.ctat-LR-fusion.fusion_predictions.tsv.gz"))

files
```


```{r}

data = NULL

for (filename in files) {
    
    message("-parsing: ", filename)
    df = fread(filename,  header=T, sep="\t", stringsAsFactors = F, drop=c("LR_accessions", "JunctionReads", "SpanningFrags", "CounterFusionLeftReads", "CounterFusionRightReads")) 
    samplename = str_replace(filename, "data/", "")
    samplename = str_replace(samplename, ".ctat-LR-fusion.fusion_predictions.tsv.gz", "")
    df$sample = samplename
    
    data = rbind(data, df)
}

data = data %>% 
    rename(fusion = `#FusionName`, SR_FFPM = FFPM) %>%  
    mutate(num_SR = est_J + est_S) %>%
    group_by(sample, fusion) %>% arrange(desc(num_LR), desc(num_SR))  %>% mutate(fusion_iso = paste(fusion, 'iso', row_number())) %>% ungroup()
  


data$lexsort_fusion_name = sapply(data$fusion, function(x) { paste(sort(str_split(x, "--")[[1]]), collapse="--") })
data = data %>% rowwise() %>% mutate(lexsort_fusion_name = paste0(sample, "|", lexsort_fusion_name))


data %>% head()

```

```{r}
# incorporate sequenced read stats:

sequenced_bases = read.table("data/DepMap.sequenced_read_stats.tsv.gz", header=T, sep="\t", stringsAsFactors = F) %>%
    select(cell_line, LR_num_bases, SR_num_bases)


sequenced_bases
```

```{r}

data = left_join(data, sequenced_bases,
                 by=c('sample'='cell_line'))

head(data)
```
```{r}

data = data %>% mutate(LR_FFPGB = num_LR / (LR_num_bases / 1e9) ,
                       SR_FFPGB = num_SR / (SR_num_bases / 1e9) )

head(data) %>% select(fusion, sample, num_LR, LR_FFPGB, num_SR, SR_FFPGB)
```


```{r}

# also, incorporate approx cost

# for $500: (very course!!)
# LR - 10e9 bases
# SR - 15e9 bases

data = data %>% mutate(LR_FFPD = num_LR / (LR_num_bases * (500/10e9)), 
                SR_FFPD = num_SR / (SR_num_bases * (500/15e9)) )



```


```{r}

write.table(data, file="DepMap_v1v2mrgd.ctatLRF_FI.consolidated.tsv", quote=F, sep="\t", row.names=F)


```




```{r}

# How many fusion gene pairs and isoforms if require both long and short read support for breakpoints?
data %>% filter(num_LR > 0 & num_SR > 0) %>% select(sample, fusion) %>% unique() %>% nrow()
# 212

data %>% filter(num_LR > 0 & num_SR > 0) %>% select(sample, fusion, LeftBreakpoint, RightBreakpoint) %>% unique() %>% nrow()
# 294

```


# compare long vs. short read fusion evidence, normalized by sequencing depth

```{r}

data %>% ggplot(aes(x=log10(LR_FFPM), y=log10(SR_FFPM))) + 
    theme_bw() +
    geom_point() +
    ggtitle("ctatLRF_FI-everything LR vs. SR fusion expression (FFPM)") +                                      
  stat_smooth(method = "lm", 
              formula = y ~ x, 
              geom = "smooth") + 
    geom_abline(intercept=0, slope=1, color='purple')


```


```{r}

cor.test(x=log2(data$LR_FFPM), y=log2(data$SR_FFPM), use='complete.obs')

# R=0.70, p<2.2e-16)

```


# compare long vs. short read fusion evidence, normalized by numbers of sequenced bases

```{r}

data %>% ggplot(aes(x=log10(LR_FFPGB), y=log10(SR_FFPGB))) + 
    theme_bw() +
    geom_point() +
    ggtitle("ctatLRF_FI-everything LR vs. SR fusion expression (FFPGB)") +                                      
  stat_smooth(method = "lm", 
              formula = y ~ x, 
              geom = "smooth") + 
    geom_abline(intercept=0, slope=1, color='purple')


```




# restrict to the TP fusions 

```{r}

TP_fusions = read.table("../3b.DepMap9Lines_Benchmarking/3b.2.Illumina_TP_unique_FP_bmarking/__illum_TP_uniq_FP.arriba,starF/data/preds.collected.gencode_mapped.wAnnot.filt.pass.proxy_assignments.byProgAgree.illum_agree.truth_set",
    #"../3b.DepMap9Lines_Benchmarking/3b.2.IncludeIlluminaSupportedFusions/data/preds.collected.gencode_mapped.wAnnot.filt.pass.proxy_assignments.byProgAgree.min_2.truth_set", 
                        header=T, sep="\t", stringsAsFactors = F) %>% 
    select(proxy_fusion_name)

TP_fusions$sample = sapply(TP_fusions$proxy_fusion_name, function(x) { str_split(x, "\\|")[[1]][1] })
 
TP_fusions$fusion = sapply(TP_fusions$proxy_fusion_name, function(x) { str_split(x, "\\|")[[1]][2] }) 

TP_fusions$lexsort_fusion_name = sapply(TP_fusions$fusion, function(x) { paste(sort(str_split(x, "--")[[1]]), collapse="--") })
TP_fusions = TP_fusions %>% rowwise() %>% mutate(lexsort_fusion_name = paste0(sample, "|", lexsort_fusion_name))

nrow(TP_fusions)

TP_fusions %>% head()
```



```{r}

data_TP_only = inner_join(TP_fusions %>% select(lexsort_fusion_name),
                  data,
                 by='lexsort_fusion_name',
                 multiple='all'
                 )



nrow(data_TP_only)

data_TP_only %>% select(sample, fusion) %>% unique() %>% nrow()

data_TP_only %>% head()

```

```{r}

# how many TP genes and isoforms have both long and short read support:
data_TP_only %>% filter(num_LR > 0 && num_SR > 0) %>% nrow()

data_TP_only %>% filter(num_LR > 0 && num_SR > 0) %>% select(sample, fusion) %>% unique() %>% nrow()


```



```{r}

# by read counts

data_TP_only %>% 
    select(sample, fusion_iso, num_LR, num_SR) %>% 
    gather(key=read_type, value=read_count, num_LR, num_SR) %>%
    ggplot(aes(x=fusion_iso, y=read_count)) + 
        theme_bw() +
        geom_bar(stat='identity', position = 'dodge', aes(fill=read_type)) +
    
        facet_grid(. ~ sample, scales = "free", space='free') +
        scale_x_discrete(expand = c(0, 0.5))  +
    
        theme(axis.text.x = element_text(angle = 90, hjust = 1))

```


```{r}

# by FFPM

data_TP_only %>% 
    select(sample, fusion_iso, LR_FFPM, SR_FFPM) %>% 
    gather(key=FFPM_type, value=FFPM, LR_FFPM, SR_FFPM) %>%
    
    mutate(FFPM = FFPM * 100) %>%
    
    ggplot(aes(x=fusion_iso, y=FFPM)) + 
        theme_bw() +
        geom_bar(stat='identity', position = 'dodge', aes(fill=FFPM_type)) +
    
    facet_grid(. ~ sample, scales = "free", space='free') +
    scale_x_discrete(expand = c(0, 0.5))  +
    
    scale_y_continuous(trans='log10') +
    
    theme(axis.text.x = element_text(angle = 90, hjust = 1)) +
    
    ylab("FFPM * 100")

```




```{r}

# label those fusions that have the most extreme difference with SR >> LR

SR_enriched_fusion_isoforms = data_TP_only %>% 
    filter(LR_FFPM > 0 & SR_FFPM > 0) %>%
    select(sample, fusion, fusion_iso, LR_FFPM, SR_FFPM, LR_FFPGB, SR_FFPGB, LR_FFPD, SR_FFPD) %>% 
    rowwise() %>%
    mutate(SR_enrichment = SR_FFPM / LR_FFPM) %>%
    mutate(SR_enrichment_GB = SR_FFPGB / LR_FFPGB) %>%
    mutate(SR_enrichment_pDollar = SR_FFPD / LR_FFPD) %>%
    arrange(desc(SR_enrichment)) %>%
    ungroup() %>% mutate(rank = row_number())


SR_enriched_fusion_isoforms %>% filter(SR_enrichment > 1)


# BCR-ABL1 has >30-fold enrichment for SR

```



```{r}

depmap_LR_vs_SR_fusion_FFPM_scatterplot = data_TP_only %>%     
    select(sample, fusion_iso, LR_FFPM, SR_FFPM) %>% 
    ggplot(aes(x=log10(LR_FFPM), y=log10(SR_FFPM))) + 
    theme_bw() +
    geom_point() +
    stat_smooth(method = "lm", 
              formula = y ~ x, 
              geom = "smooth") + 
    geom_abline(intercept=0, slope=1, color='purple') +
    geom_point(data=SR_enriched_fusion_isoforms %>% filter(SR_enrichment>=3), color='red')  +
    geom_text_repel(data=SR_enriched_fusion_isoforms %>% filter(SR_enrichment>=3), aes(label=fusion_iso)) +
    ggtitle("CTAT-LR-FI FFPM Comparison for isoforms of TP fusions")


depmap_LR_vs_SR_fusion_FFPM_scatterplot

```

```{r}

ggsave(depmap_LR_vs_SR_fusion_FFPM_scatterplot, file="depmap_LR_vs_SR_fusion_FFPM_scatterplot.svg", width=6, height=5)


```


```{r}


cor.test(x=log2(data_TP_only$LR_FFPM), y=log2(data_TP_only$SR_FFPM), use='complete.obs')

# with TP fusions:  R=0.66, p<2.2e-16

```


```{r}

# examine according to num sequenced bases:

data_TP_only %>%     
    select(sample, fusion_iso, LR_FFPGB, SR_FFPGB) %>% 
    ggplot(aes(x=log10(LR_FFPGB), y=log10(SR_FFPGB))) + 
    theme_bw() +
    geom_point() +
    stat_smooth(method = "lm", 
              formula = y ~ x, 
              geom = "smooth") + 
    geom_abline(intercept=0, slope=1, color='purple') +
    geom_point(data=SR_enriched_fusion_isoforms %>% filter(SR_enrichment>=3), color='red')  +
    geom_text_repel(data=SR_enriched_fusion_isoforms %>% filter(SR_enrichment>=3), aes(label=fusion_iso)) +
    ggtitle("CTAT-LR-FI FFPGB Comparison for isoforms of TP fusions")



```


```{r}

# examine enrichment by cost

data_TP_only %>%     
    select(sample, fusion_iso, LR_FFPD, SR_FFPD) %>% 
    ggplot(aes(x=log10(LR_FFPD), y=log10(SR_FFPD))) + 
    theme_bw() +
    geom_point() +
    stat_smooth(method = "lm", 
              formula = y ~ x, 
              geom = "smooth") + 
    geom_abline(intercept=0, slope=1, color='purple') +
    geom_point(data=SR_enriched_fusion_isoforms %>% filter(SR_enrichment>=3), color='red')  +
    geom_text_repel(data=SR_enriched_fusion_isoforms %>% filter(SR_enrichment>=3), aes(label=fusion_iso)) +
    ggtitle("CTAT-LR-FI FFpDollar Comparison for isoforms of TP fusions")





```




```{r}

depmap_SR_enrichment_rankings_plot = SR_enriched_fusion_isoforms %>%
    mutate(rn = row_number() ) %>%
    ggplot(aes(x=rn, y=SR_enrichment)) + 
    theme_bw() +
    geom_point() + geom_hline(yintercept = 1.0, color='purple') +
    scale_y_continuous(trans='log10') +
    xlab("Fusion isoform ranked by SR_enrichment") +
    ggtitle("SR/LR enrichment by FFPM")
   
depmap_SR_enrichment_rankings_plot

```


```{r}

#ggsave(depmap_SR_enrichment_rankings_plot, file="depmap_SR_enrichment_rankings_plot.svg", width=5, height=3)


```


```{r}
# redo using normalization by sequencing depth:


SR_enriched_fusion_isoforms %>% 
    arrange(desc(SR_enrichment_GB)) %>%
    mutate(rn = row_number() ) %>%
    ggplot(aes(x=rn, y=SR_enrichment_GB)) + 
    theme_bw() +
    geom_point() + geom_hline(yintercept = 1.0, color='purple') +
    scale_y_continuous(trans='log10') +
    xlab("Fusion isoform ranked by SR_enrichment_GB") +
    ggtitle("SR/LR enrichment by FFPGB")


```


```{r}

median_SR_enrichment_GB = SR_enriched_fusion_isoforms %>% summarize(median(SR_enrichment_GB))

median_SR_enrichment_GB 
```

```{r}
# redo using normalization by sequencing depth:

SR_enriched_fusion_isoforms = SR_enriched_fusion_isoforms %>% 
    arrange(desc(SR_enrichment_GB)) %>%
    mutate(rn = row_number() ) 


SR_vs_LR_FFPGB_ranking_plot = SR_enriched_fusion_isoforms %>%
    ggplot(aes(x=rn, y=SR_enrichment_GB)) + 
    theme_bw() +
    geom_point() + 
    geom_hline(yintercept = 1.0, color='purple') +
    geom_hline(yintercept = median_SR_enrichment_GB$`median(SR_enrichment_GB)`, color='blue')  +
    scale_y_continuous(trans='log10') +
    geom_point(data=SR_enriched_fusion_isoforms %>% filter(SR_enrichment>=3), color='red')  +
    geom_text_repel(data=SR_enriched_fusion_isoforms %>% filter(SR_enrichment>=3), aes(label=fusion_iso)) +
    xlab("Fusion isoform ranked by SR_enrichment_GB") +
    ggtitle("SR/LR enrichment by FFPGB")


SR_vs_LR_FFPGB_ranking_plot

```

```{r}

ggsave(SR_vs_LR_FFPGB_ranking_plot, filename="SR_vs_LR_FFPGB_ranking_plot.svg", width=6, height=5)

```




```{r}
# examine by approximate cost.

SR_enriched_fusion_isoforms %>% 
    arrange(desc(SR_enrichment_pDollar)) %>%
    mutate(rn = row_number() ) %>%
    ggplot(aes(x=rn, y=SR_enrichment_pDollar)) + 
    theme_bw() +
    geom_point() + geom_hline(yintercept = 1.0, color='purple') +
    scale_y_continuous(trans='log10') +
    xlab("Fusion isoform ranked by SR_enrichment_per$") +
    ggtitle("SR/LR enrichment by FFPD")


```



# examine specific fusions and breakpoint splice support

```{r}

plot_fusion_expression_by_breakpoint = function(sample_name, fusion_name) {
    
    df = data %>% filter(sample == sample_name & fusion == fusion_name) %>% 
        select(sample, fusion, LeftLocalBreakpoint, RightLocalBreakpoint, num_LR, est_J, LR_FFPM, SR_FFPM) %>% 
        arrange(sample, RightLocalBreakpoint, LeftLocalBreakpoint) %>%
        mutate(num_LR = ifelse (is.na(num_LR), 0, num_LR)) %>%
        mutate(LR_FFPM = ifelse(is.na(LR_FFPM), 0, LR_FFPM)) %>%
        mutate(est_J = ifelse(is.na(est_J), 0, est_J)) 
    
    print(df)
    
    
    p1 = df %>% gather(key=readtype, value=readcount, num_LR, est_J) %>% 
        filter(readcount > 0) %>%
        ggplot(aes(x=RightLocalBreakpoint, y=LeftLocalBreakpoint) ) + 
        theme_bw() +
        ggtitle(paste(sample_name, fusion_name, "by read count")) +
        geom_point(aes(color=readtype, size=readcount), alpha=0.5) +
        facet_wrap(~readtype) 
    
    plot(p1)
    
    
    p2 = df %>% gather(key=readtype, value=FFPM, LR_FFPM, SR_FFPM) %>% 
        filter(FFPM > 0) %>%
        ggplot(aes(x=RightLocalBreakpoint, y=LeftLocalBreakpoint) ) + 
        theme_bw() +
        ggtitle(paste(sample_name, fusion_name, "by FFPM")) +
        geom_point(aes(color=readtype, size=FFPM), alpha=0.5) +
        facet_wrap(~readtype) 
    
    plot(p2)
    
    p3 = df %>% gather(key=readtype, value=FFPM, LR_FFPM, SR_FFPM) %>% 
        ggplot(aes(x=RightLocalBreakpoint, y=LeftLocalBreakpoint) ) + 
        theme_bw() +
        ggtitle(paste(sample_name, fusion_name, "by FFPM")) +
        geom_point(aes(color=readtype, fill=readtype, size=FFPM), alpha=0.5) 
    
    plot(p3)
    
    
    p4 = df %>% gather(key=readtype, value=FFPM, LR_FFPM, SR_FFPM) %>% 
        mutate(breakpoint = paste(RightLocalBreakpoint, LeftLocalBreakpoint)) %>% 
        ggplot(aes(x=breakpoint, y=FFPM, fill=readtype) ) + 
        theme_bw() +
        ggtitle(paste(sample_name, fusion_name, "by FFPM")) +
        geom_bar(stat='identity', position='dodge')  +
        theme(axis.text.x = element_text(angle = 90, hjust = 1))
    
    plot(p4)
    
    lm_df = lm(df$SR_FFPM ~ df$LR_FFPM)
    print(lm_df)
    lm_y_intercept=lm_df$coefficients[1]
    lm_slope=lm_df$coefficients[2]
    
    
     p5 = df %>%
        mutate(breakpoint = paste(RightLocalBreakpoint, LeftLocalBreakpoint)) %>% 
        ggplot(aes(x=LR_FFPM, y=SR_FFPM) ) +  
        theme_bw() + 
        geom_point() +  geom_abline(slope=lm_slope, intercept = lm_y_intercept) +
        ggtitle(paste(sample_name, fusion_name, "by FFPM")) 
    
    plot(p5)
    
    if(nrow(df)>2) {
        print(cor.test(df$LR_FFPM, df$SR_FFPM))
    }
}

```




## TATDN1--GSDMB
```{r}

sample_name = "SKBR3"
fusion_name = "TATDN1--GSDMB"

plot_fusion_expression_by_breakpoint(sample_name, fusion_name)

```

## K562 BCR--ABL1

```{r}



plot_fusion_expression_by_breakpoint("K562", "BCR--ABL1")


```



```{r}

# HCC1187	PUM1--TRERF1 

plot_fusion_expression_by_breakpoint("HCC1187", "PUM1--TRERF1")


```



```{r}

# VCAP	TMPRSS2--ERG

plot_fusion_expression_by_breakpoint("VCAP", "TMPRSS2--ERG")


```


# 
```{r}

plot_fusion_expression_by_breakpoint("SKBR3", "SAMD12--MRPL13")

```




# Examine all fusions with multiple isoforms

How many fusions do we have multiple isoforms with evidence?  How many with both short and long read support?

Are expression values correlated among isoforms according to read types?

Is the dominant isoform the same?



```{r}

# how many total TP fusions did ctatLRF find evidence for with both short and long reads?

fusions_w_both_long_and_short = data %>% 
    filter(num_LR > 0 & num_SR > 0) %>%
   select(sample, fusion) %>% unique()

nrow(fusions_w_both_long_and_short)

```

```{r}

min_multi_isoforms = 3

```



```{r}

# how many have multiple fusion isoforms?

mult_isoform_data = left_join(fusions_w_both_long_and_short, data, 
                              by=c('sample', 'fusion'),
                              multiple='all') %>% 
    group_by(sample, fusion) %>% filter(n()>=min_multi_isoforms) %>% ungroup()

mult_isoform_data %>% select(sample, fusion) %>% unique() %>% nrow()


```



```{r}

# how many have both LR and SR isoform support?

mult_isoform_data_both_read_types = mult_isoform_data %>% filter(num_LR > 0 & est_J > 0) %>% group_by(sample, fusion) %>% filter(n()>min_multi_isoforms) %>% ungroup() 

mult_iso_both_sample_fusion_names = mult_isoform_data_both_read_types %>% group_by(sample, fusion) %>% tally(name='num_multi_iso_both_types')

# include isoformws supported uniquely by long or short in the downstream analysis here:
mult_isoform_data_both_read_types = left_join(mult_iso_both_sample_fusion_names, mult_isoform_data, by=c('sample', 'fusion'))

mult_iso_both_sample_fusion_names %>% arrange(desc(num_multi_iso_both_types))

# 6 fusions with at least 3 splicing isoforms each:

```


## Compare read support for fusion isoforms across each fusion gene


```{r}


depmap_fusion_isoform_expr_LR_SR_comparison_plot =  mult_isoform_data_both_read_types %>% 
    ggplot(aes(x=LR_FFPM, y=SR_FFPM)) + 
    theme_bw() +
    geom_point(aes(color=fusion)) +
    ggtitle("restricted to fusion genes w/ multi isoforms supported by both read types") +
    stat_smooth(method = "lm", 
              formula = y ~ x, 
              geom = "smooth") +
    facet_wrap(~fusion, scale='free')

depmap_fusion_isoform_expr_LR_SR_comparison_plot
```


```{r}


depmap_fusion_isoform_expr_LR_SR_comparison_plot =  mult_isoform_data_both_read_types %>% 
    ggplot(aes(x=log10(LR_FFPM), y=log10(SR_FFPM))) + 
    theme_bw() +
    geom_point(aes(color=fusion)) +
    ggtitle("restricted to fusion genes w/ multi isoforms supported by both read types") +
    stat_smooth(method = "lm", 
              formula = y ~ x, 
              geom = "smooth") +
    facet_wrap(~fusion, scale='free')

depmap_fusion_isoform_expr_LR_SR_comparison_plot
```

```{r}

ggsave(depmap_fusion_isoform_expr_LR_SR_comparison_plot, file="depmap_fusion_isoform_expr_LR_SR_comparison_plot.svg", width=9, height=5)


```


```{r}

for (i in seq(nrow(mult_iso_both_sample_fusion_names))) {
    
    sample_name = mult_iso_both_sample_fusion_names[i,1]$sample
    fusion_name =  mult_iso_both_sample_fusion_names[i,2]$fusion
    
    loc_df = mult_isoform_data_both_read_types %>% 
        filter(sample==sample_name & fusion==fusion_name) %>%
        select(sample, fusion, LeftLocalBreakpoint, RightLocalBreakpoint, num_LR, est_J, LR_FFPM, SR_FFPM) %>% 
        arrange(sample, RightLocalBreakpoint, LeftLocalBreakpoint) 
    
    #print(loc_df)
    
    lm_loc_df = lm(loc_df$SR_FFPM ~ loc_df$LR_FFPM)
    print(lm_loc_df)
    lm_y_intercept=lm_loc_df$coefficients[1]
    lm_slope=lm_loc_df$coefficients[2]
    
    p = loc_df %>%
        mutate(breakpoint = paste(RightLocalBreakpoint, LeftLocalBreakpoint)) %>% 
        ggplot(aes(x=LR_FFPM, y=SR_FFPM) ) + 
        theme_bw() +
        geom_point() + geom_abline(slope=lm_slope, intercept = lm_y_intercept) +
        ggtitle(paste(sample_name, fusion_name, "by FFPM")) 
    
    plot(p)
    
    if(nrow(loc_df)>2) {
        message(fusion_name, " Cor: FFPM  ")
        print(cor.test(loc_df$LR_FFPM, loc_df$SR_FFPM))
    
        message(fusion_name, " Cor: log10(FFPM)  ")
        print(cor.test(log10(loc_df$LR_FFPM), log10(loc_df$SR_FFPM)))
    }
    
    
}


```

```{r}

plot_fusion_expression_by_breakpoint("SKBR3", "CYTH1--EIF3H")

# significant correlation:  R=0.997, p=1.9e-4

```



# include gene structure info in plots

```{r}

fusion_gene_structures = read.table("data/FI.exons.dat", header=F, sep="\t")
colnames(fusion_gene_structures) = c('fusion_name', 'feature_type', 'lend', 'rend')

```



```{r}


plot_fusion_expression_by_breakpoint_incl_gene_structures = function(sample_name, fusion_name) {
    
    
    ret_plots = list()
    
    
    df = data %>% filter(sample == sample_name & fusion == fusion_name) %>% 
        select(sample, fusion, LeftLocalBreakpoint, RightLocalBreakpoint, num_LR, est_J, LR_FFPM, SR_FFPM) %>% 
        arrange(sample, RightLocalBreakpoint, LeftLocalBreakpoint) %>%
        mutate(num_LR = ifelse (is.na(num_LR), 0, num_LR)) %>%
        mutate(LR_FFPM = ifelse(is.na(LR_FFPM), 0, LR_FFPM)) %>%
        mutate(est_J = ifelse(is.na(est_J), 0, est_J)) 
    
    print(df)
    
    
    p = df %>% gather(key=readtype, value=FFPM, LR_FFPM, SR_FFPM) %>% 
        filter(FFPM > 0) %>%
        ggplot() +
        theme_bw() +
        ggtitle(paste(sample_name, fusion_name, "by FFPM")) +
        geom_point(aes(x=LeftLocalBreakpoint, y=RightLocalBreakpoint, color=readtype, size=FFPM), alpha=0.5)
    #+
    #    facet_wrap(~readtype) 
    
    fname = fusion_name
    
    gene_structure_info = fusion_gene_structures %>% filter(fusion_name == fname)
    
    geneA_info = gene_structure_info %>% filter(feature_type == "GeneA")
    geneB_info = gene_structure_info %>% filter(feature_type == "GeneB")
    
    print(geneA_info)
    print(geneB_info)

    ## plot gene structures.
    padding = 1000

    geneA_lend = min(geneA_info$lend)
    geneA_rend = max(geneA_info$rend)
    geneA_length = geneA_rend - geneA_lend + 2*padding

    geneB_lend = min(geneB_info$lend)
    geneB_rend = max(geneB_info$rend)
    geneB_length = geneB_rend - geneB_lend + 2*padding

    

    ## draw geneA
    geneA_minY = geneB_lend - 0.95*padding
    geneA_maxY = geneA_minY + 0.03*geneB_length

    p1 = p + geom_rect(data=geneA_info, aes(xmin=lend, xmax=rend, ymin=geneA_minY, ymax=geneA_maxY), fill=NA, color='black', alpha=0.5)

    geneA_midY = mean(c(geneA_minY, geneA_maxY))

    p1 = p1 + geom_segment(data=geneA_info, aes(x=geneA_lend, xend=geneA_rend, y=geneA_midY, yend=geneA_midY), alpha=0.5) # geneB center line

    ## draw geneB

    geneB_minX = geneA_lend - 0.95*padding
    geneB_maxX = geneB_minX + 0.03*geneA_length

    p1 = p1 + geom_rect(data=geneB_info, aes(ymin=lend, ymax=rend, xmin=geneB_minX, xmax=geneB_maxX), fill=NA, color='black', alpha=0.5)

    geneB_midX = mean(c(geneB_minX, geneB_maxX))
    p1 = p1 + geom_segment(data=geneB_info, aes(x=geneB_midX, xend=geneB_midX, y=geneB_lend, yend=geneB_rend), alpha=0.5) # geneB center line


    
    ## set up scales for view
    p1 = p1 +
        xlim(geneA_lend - padding, geneA_rend + padding) +
        ylim(geneB_lend - padding, geneB_rend + padding)

    
    plot(p1)
    ret_plots[[1]] = p1
    
    p1 = p1 + facet_wrap(~readtype) 
    plot(p1)
    ret_plots[[2]] = p1
    
    
    # zoom in around breakpoints
    
    min_right_breakpoint = min(df$RightLocalBreakpoint)
    max_right_breakpoint = max(df$RightLocalBreakpoint)
    
    min_left_breakpoint = min(df$LeftLocalBreakpoint)
    max_left_breakpoint = max(df$LeftLocalBreakpoint)
    
    
    padding = 0.25 * padding
    
     ## draw geneA (along X-axis)
    geneA_minY = min_right_breakpoint - 0.95*padding
    geneA_maxY = geneA_minY + 0.03* (max_right_breakpoint - min_right_breakpoint + 2*padding)

    p2 = p + geom_rect(data=geneA_info, aes(xmin=lend, xmax=rend, ymin=geneA_minY, ymax=geneA_maxY), fill=NA, color='black', alpha=0.5)

    geneA_midY = mean(c(geneA_minY, geneA_maxY))

    p2 = p2 + geom_segment(data=geneA_info, aes(x=min_left_breakpoint-padding, xend=max_left_breakpoint+padding, y=geneA_midY, yend=geneA_midY), alpha=0.5) # geneB center line

    ## draw geneB (along Y-axis)

    geneB_minX = min_left_breakpoint - 0.95*padding
    geneB_maxX = geneB_minX + 0.03* (max_left_breakpoint - min_left_breakpoint + 2*padding)

    p2 = p2 + geom_rect(data=geneB_info, aes(ymin=lend, ymax=rend, xmin=geneB_minX, xmax=geneB_maxX), fill=NA, color='black', alpha=0.5)

    geneB_midX = mean(c(geneB_minX, geneB_maxX))
    p2 = p2 + geom_segment(data=geneB_info, aes(x=geneB_midX, xend=geneB_midX, y=min_right_breakpoint-padding, yend=max_right_breakpoint+padding), alpha=0.5) # geneB center line

    
    p2 = p2 + ylim(min_right_breakpoint-padding, max_right_breakpoint+padding) +
        xlim(min_left_breakpoint-padding, max_left_breakpoint+padding)
    
    plot(p2)
    ret_plots[[3]] = p2
    
    p2 = p2 + facet_wrap(~readtype)
    plot(p2)
    ret_plots[[4]] = p2
    
    
    return(ret_plots)
    
}

```


```{r}

plot_fusion_expression_by_breakpoint_incl_gene_structures("K562", "BCR--ABL1")

```


```{r}

plots = plot_fusion_expression_by_breakpoint_incl_gene_structures("VCAP", "TMPRSS2--ERG")

ggsave(plots[[2]], file="depmap_VCaP_TMPRSS2--ERG_LR_vs_SR_isoforms.svg", width=9, height=6)

ggsave(plots[[4]], file="depmap_VCaP_TMPRSS2--ERG_LR_vs_SR_isoforms.zoomed.svg", width=9, height=6)

# has 1 isoform uniquely supported by long reads

```


```{r}


plots = plot_fusion_expression_by_breakpoint_incl_gene_structures("SKBR3", "CYTH1--EIF3H")

ggsave(plots[[1]], file="depmap_SKBR3_CYTH1--EIF3H_LR_vs_SR_isoforms.svg", width=9, height=6)


ggsave(plots[[3]], file="depmap_SKBR3_CYTH1--EIF3H_LR_vs_SR_isoforms.zoomed.svg", width=6, height=4)

# save the structure fore the full-length view of each gene:
ggsave(plots[[1]] + xlim(800,1000) + theme_bw(), file="EIF3H_structure.svg", width=4)


```


```{r}


cyth1_eif3h_cor_plot =  mult_isoform_data_both_read_types %>% 
    filter(fusion=="CYTH1--EIF3H") %>%
    ggplot(aes(x=log10(LR_FFPM), y=log10(SR_FFPM))) + 
    theme_bw() +
    geom_point(aes(color=fusion), size=rel(3), color='black') +
    ggtitle("restricted to fusion genes w/ multi isoforms supported by both read types") +
    stat_smooth(method = "lm", 
              formula = y ~ x, 
              geom = "smooth") 
cyth1_eif3h_cor_plot 

ggsave(cyth1_eif3h_cor_plot, file="cyth1_eif3h_LR_SR_correlation_plot.svg", width=5, height=4)
```



```{r}


plots = plot_fusion_expression_by_breakpoint_incl_gene_structures("SKBR3", "TATDN1--GSDMB")

ggsave(plots[[2]], file="depmap_SKBR3_TATDN1--GSDMB_LR_vs_SR_isoforms.svg", width=9, height=6)

ggsave(plots[[4]], file="depmap_SKBR3_TATDN1--GSDMB_LR_vs_SR_isoforms.zoomed.svg", width=9, height=6)

# has 4 isoforms uniquely supported by SR

```









