---
title: "CTAT_SeraCareFusion"
author: "bhaas"
date: '2024-02-03'
output: github_document
---

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


Fusion prediction results and benchmarking for SeraCare fusion v4 mix are available here:
https://github.com/fusiontranscripts/LR-FusionBenchmarking/tree/master/SeraCareFusions

Analysis performed below:

The Iso-seq labeled sample is actually MAS-Iso-seq 'monomer' type.
The other two MAS-Iso-seq used the 8-concatemer arrays followed by deconcatenation.


# number of fusion reads for each control fusion

```{r}

control_fusions = read.table("data/SeraCare_fusion_targets.tsv", header=T, sep="\t")

control_fusions 
```


```{r}

# read in the ctat-LRF fusion predictions for each of the three samples

rep1_data = read.table("data/ctatLRF_FI/SeraCarev4Fusion.IsoSeq.m64020e_230609_061545.refined.Oct32023.ctat-LR-fusion.fusion_predictions.tsv.gz", header=T, com='', sep="\t") %>%
    mutate(dataset= 'ISO-seq')

rep2_data = read.table("data/ctatLRF_FI/SeraCarev4Fusion.m64363e_230614_195212.refine.bc08.MASseq.ctat-LR-fusion.fusion_predictions.tsv.gz", header=T, com='', sep="\t") %>%
    mutate(dataset = 'MAS-seq-R1')

rep3_data = read.table("data/ctatLRF_FI/SeraCarev4Fusion.m64363e_230614_195212.refine.bc09.MASseq.ctat-LR-fusion.fusion_predictions.tsv.gz", header=T, com='', sep="\t") %>%
    mutate(dataset = 'MAS-seq-R2')


ctatLRF_FI_data = bind_rows(rep1_data, rep2_data, rep3_data) %>% rename(FusionName = X.FusionName, SR_FFPM = FFPM) %>%
    mutate(num_SR = est_J + est_S) 

ctatLRF_FI_data %>% head()


```

```{r}

p = ctatLRF_FI_data %>% ggplot(aes(x=LR_FFPM, y=SR_FFPM)) + geom_point() +  
    theme_bw() +
 stat_smooth(method = "lm", 
              formula = y ~ x, 
              geom = "smooth")

p
```

```{r}

lm(ctatLRF_FI_data$LR_FFPM ~ ctatLRF_FI_data$SR_FFPM  )

```


```{r}

p + geom_point(data=ctatLRF_FI_data %>% filter(FusionName %in% control_fusions$FusionName), color='red')

```


```{r}
# conpare FFPM values LR vs. SR

ctatLRF_FI_data %>% select(FusionName, LR_FFPM, SR_FFPM) %>% 
    filter (! is.na(SR_FFPM)) %>% filter(! is.na(LR_FFPM)) %>%
    mutate(LR_div_SR_FFPM = LR_FFPM/SR_FFPM) %>% 
    mutate(control_fusion = (FusionName %in% control_fusions$FusionName)) %>%
    ggplot() + 
    theme_bw() +
    geom_boxplot(aes(x=control_fusion, y=LR_div_SR_FFPM))


```

The spike-in control fusions do appear to be better captured by MAS-ISO-seq than the non-control fusions.



```{r}

# restrict fusion calls to the 16 control fusions

ctatLRF_FI_control_results = left_join( cross_join(control_fusions, data.frame(dataset=c('ISO-seq', 'MAS-seq-R1', 'MAS-seq-R2'))),
                                    ctatLRF_FI_data %>% 
                                        select(FusionName, LeftBreakpoint, RightBreakpoint, 
                                               num_LR, LR_FFPM, num_SR, SR_FFPM, dataset),
                                    by=c('FusionName', 'LeftBreakpoint', 'RightBreakpoint', 'dataset'))


ctatLRF_FI_control_results =  ctatLRF_FI_control_results %>% mutate(num_SR = ifelse(is.na(num_SR), 0, num_SR))

ctatLRF_FI_control_results

```

# examine number of long reads per control fusion

```{r}

# by read count

seracare_read_count_barplot = ctatLRF_FI_control_results %>% 
    ggplot(aes(x=FusionName, y=num_LR)) + theme_bw() +
    geom_bar(stat='identity', position='dodge', aes(fill=dataset))   + 
    theme(axis.text.x = element_text(angle = 90, hjust = 1))

seracare_read_count_barplot

```

```{r}

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

```

```{r}

# by normalized read count

seracare_FFPM_barplot = ctatLRF_FI_control_results %>% 
    ggplot(aes(x=FusionName, y=LR_FFPM)) + theme_bw() +
    geom_bar(stat='identity', position='dodge', aes(fill=dataset))   + 
    theme(axis.text.x = element_text(angle = 90, hjust = 1))


seracare_FFPM_barplot

```

```{r}

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

```

```{r}

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

```

# examine counts of long reads (ctatLRF) vs. short reads (FI) evidence

```{r}

paperfig_LR_vs_SR = ctatLRF_FI_control_results %>% gather(key=read_count_type, value=read_count, num_LR, num_SR) %>%
    ggplot(aes(x=read_count_type, y=read_count)) +  theme_bw() +
    geom_bar(stat='identity', position = 'dodge', aes(color=dataset, fill=read_count_type)) +
    facet_wrap(~FusionName) +
      theme(axis.text.x = element_text(angle = 90, hjust = 1))  +
          geom_text(position=position_dodge(width=0.9), aes(label=read_count, group=dataset))

paperfig_LR_vs_SR

```


```{r}

ggsave(paperfig_LR_vs_SR, file="seracare_LR_vs_SR_barplot.svg", width=9, height=7)

```


```{r}
# How many fusions are missing short read evidence in at least one of the replicates?

ctatLRF_FI_control_results %>% filter(num_SR == 0) %>% select(FusionName) %>% unique()

# for 9/16 control fusions, missing SR evidence in at least one replicate

```

```{r}
# what's the range of long read support observed?

ctatLRF_FI_control_results %>% summarize(max_LR_support=max(num_LR), min_LR_support = min(num_LR))

# 20 to 139 long read support for control fusions
```



```{r}

# compare the normalized read support FFPM values

ctatLRF_FI_control_results %>% ggplot(aes(x=LR_FFPM, y=SR_FFPM)) + geom_point() +
    theme_bw() +
  stat_smooth(method = "lm", 
              formula = y ~ x, 
              geom = "smooth")


```


```{r}


lm(ctatLRF_FI_control_results$SR_FFPM ~ ctatLRF_FI_control_results$LR_FFPM)

```

```{r}

plot(ctatLRF_FI_control_results$LR_FFPM, ctatLRF_FI_control_results$SR_FFPM )
abline(a=0.05671, b=0.02847)


```

Based on the slope of the line, would be ~35x better detection based on LR than SR, but obviously a lot of variance around that.



```{r}

stats::fivenum(ctatLRF_FI_control_results$LR_FFPM)

```

```{r}

stats::fivenum(ctatLRF_FI_control_results$SR_FFPM)

```

```{r}

stats::fivenum(ctatLRF_FI_control_results$LR_FFPM) / stats::fivenum(ctatLRF_FI_control_results$SR_FFPM)

```

so ~20x to ~35x more LR than SR for control fusion reads.


```{r}
# what about the non-control fusions?

noncontrol_fusions = ctatLRF_FI_data %>% filter(! FusionName %in% control_fusions$FusionName) %>% filter(LR_FFPM > 0 & SR_FFPM > 0)

noncontrol_fusions 
```

```{r}


noncontrol_fusions$LR_FFPM

noncontrol_fusions$SR_FFPM

noncontrol_fusions$LR_FFPM / noncontrol_fusions$SR_FFPM

stats::fivenum(noncontrol_fusions$LR_FFPM / noncontrol_fusions$SR_FFPM)
```

Don't have many good examples, but rougly 1x to 10x LR to SR for LR support seems likely 




# compare gene expression levels between long and short reads:

```{r}

LR_gene_expr = read.table("data/expr_quant/LR.salmon.gene.TMM.EXPR.matrix.gz", header=T, sep="\t") %>%
    rename(gene_id=X)

TruSeq_gene_expr = read.table("data/expr_quant/TruSeq.salmon.gene.TMM.EXPR.matrix.gz", header=T, sep="\t") %>%
    rename(gene_id=X)

LR_mean_gene_expr = LR_gene_expr %>% rowwise() %>%  mutate(mean_LR_gene_expr = mean(SCF_IsoSeq, SCF_bc08, SCF_bc09)) %>% select(gene_id, mean_LR_gene_expr)

TruSeq_mean_gene_expr = TruSeq_gene_expr %>% rowwise() %>% mutate(mean_TruSeq_gene_expr = mean(SeraFusion_r1, SeraFusion_r2, SeraFusion_r3)) %>% select(gene_id, mean_TruSeq_gene_expr)

mean_gene_expr = full_join(LR_mean_gene_expr, TruSeq_mean_gene_expr, by='gene_id') %>% filter(mean_LR_gene_expr > 0 & mean_TruSeq_gene_expr > 0)

```


```{r}

# from: https://stackoverflow.com/questions/13094827/how-to-reproduce-smoothscatters-outlier-plotting-in-ggplot


densVals <- function(x, y = NULL, nbin = 128, bandwidth, range.x) {
  dat <- cbind(x, y)
  # limit dat to strictly finite values
  sel <- is.finite(x) & is.finite(y)
  dat.sel <- dat[sel, ]
  # density map with arbitrary graining along x and y
  map   <- grDevices:::.smoothScatterCalcDensity(dat.sel, nbin, bandwidth)
  map.x <- findInterval(dat.sel[, 1], map$x1)
  map.y <- findInterval(dat.sel[, 2], map$x2)
  # weighted mean of the fitted density map according to how close x and y are
  # to the arbitrary grain of the map
  den <- mapply(function(x, y) weighted.mean(x = c(
    map$fhat[x, y], map$fhat[x + 1, y + 1],
    map$fhat[x + 1, y], map$fhat[x, y + 1]), w = 1 / c(
    map$x1[x] + map$x2[y], map$x1[x + 1] + map$x2[y + 1],
    map$x1[x + 1] + map$x2[y], map$x1[x] + map$x2[y + 1])),
    map.x, map.y)
  # replace missing density estimates with NaN
  res <- rep(NaN, length(sel))
  res[sel] <- den
  res
}

mean_gene_expr$point_density <- densVals(log10(mean_gene_expr$mean_LR_gene_expr), log10(mean_gene_expr$mean_TruSeq_gene_expr))

mean_gene_expr_plot = mean_gene_expr  %>% 
    ggplot(aes(x = log10(mean_LR_gene_expr), y=log10(mean_TruSeq_gene_expr))) +
    theme_bw() +
  stat_density2d(geom = "raster", aes(fill = ..density.. ^ 0.25), contour = FALSE) +
  #scale_x_log10() + scale_y_log10() +
  scale_fill_gradientn(colours = colorRampPalette(c("white", blues9))(256)) +
  # select only the 100 sparesest points
   geom_point(data = mean_gene_expr %>% arrange(point_density) %>% head(100), size = .5) +
  stat_smooth(method = "lm", 
              formula = y ~ x, 
              geom = "smooth")


mean_gene_expr_plot

```

```{r}

ggsave(mean_gene_expr_plot, file="SeraCare_LR_vs_SR_gene_expr_scatter_density_plot.svg", width=7, height=4)



```



```{r}


cor.test(log10(mean_gene_expr$mean_LR_gene_expr), log10(mean_gene_expr$mean_TruSeq_gene_expr))


```

# examine Arriba short read findings:

```{r}

arriba_1 = read.csv("data/arriba/MDL_163_SeraSeqFusion_1.arriba.fusions.tsv.gz", header=T, sep="\t", com='', stringsAsFactors = F)
arriba_1$dataset = "ISO-seq"


arriba_2 = read.csv("data/arriba/MDL_163_SeraSeqFusion_2.arriba.fusions.tsv.gz", header=T, sep="\t", com='', stringsAsFactors = F)
arriba_2$dataset = "MAS-seq-R1"


arriba_3 = read.csv("data/arriba/MDL_163_SeraSeqFusion_2.arriba.fusions.tsv.gz", header=T, sep="\t", com='', stringsAsFactors = F)
arriba_3$dataset = "MAS-seq-R2"


arriba_data = bind_rows(arriba_1, arriba_2, arriba_3)

arriba_data = arriba_data %>% 
    mutate(FusionName = paste0(X.gene1, "--", gene2) ) %>%
    mutate(num_reads = split_reads1 + split_reads2 + discordant_mates) %>%
    select(FusionName, num_reads, dataset)

arriba_data %>% head()
```

```{r}

ctatLRF_FI_control_results_w_arriba = left_join(ctatLRF_FI_control_results  %>% filter(FusionName %in% control_fusions$FusionName),
                                                arriba_data %>% filter(FusionName %in% control_fusions$FusionName) %>% mutate(arriba_SR = num_reads),
                                                by=c('FusionName', 'dataset') ) %>%
    mutate(arriba_SR = ifelse(is.na(arriba_SR), 0, arriba_SR))

```



```{r}

ctatLRF_FI_control_results_w_arriba %>% gather(key=read_count_type, value=read_count, num_LR, num_SR, arriba_SR) %>%
    mutate(read_count_type = factor(read_count_type, levels=c('num_LR', 'num_SR', 'arriba_SR'))) %>%
    ggplot(aes(x=read_count_type, y=read_count)) +  theme_bw() +
    geom_bar(stat='identity', position = 'dodge', aes(color=dataset, fill=read_count_type)) +
    facet_wrap(~FusionName) +
      theme(axis.text.x = element_text(angle = 90, hjust = 1))  +
          geom_text(position=position_dodge(width=0.9), aes(label=read_count, group=dataset))

```


# Examine reads per gene in LR vs. SR for all genes

```{r}

LR_read_gene_counts = read.csv("data/expr_quant/LR.salmon.gene.counts.matrix.gz", header=T, sep="\t", stringsAsFactors = F) %>%
    rename(gene_id = X)

colnames(LR_read_gene_counts) = c('gene_id', 'LR_r1', 'LR_r2', 'LR_r3')

LR_read_gene_counts %>% head()


```


```{r}

SR_read_gene_counts = read.csv("data/expr_quant/TruSeq.salmon.gene.counts.matrix.gz", header=T, sep="\t", stringsAsFactors = F) %>%
    rename(gene_id = X)

colnames(SR_read_gene_counts) = c('gene_id', 'SR_r1', 'SR_r2', 'SR_r3')

SR_read_gene_counts %>% head()
```


```{r}

gene_read_counts = left_join(LR_read_gene_counts, SR_read_gene_counts, by='gene_id')

gene_read_counts %>% head()
```

```{r}

gene_read_counts %>% ggplot(aes(x=LR_r1, y=SR_r1)) + geom_point() +
    theme_bw() + 
    scale_x_continuous(trans='log10') +
    scale_y_continuous(trans='log10') +
  stat_smooth(method = "lm", 
              formula = y ~ x, 
              geom = "smooth")

```

```{r}


# convert to CPM

gene_ids = gene_read_counts %>% pull(gene_id)

col_sums = gene_read_counts %>% select(-gene_id) %>% summarize(sum_LR_r1 = sum(LR_r1), 
                                                               sum_LR_r2 = sum(LR_r2),
                                                               sum_LR_r3 = sum(LR_r3),
                                                               sum_SR_r1 = sum(SR_r1),
                                                               sum_SR_r2 = sum(SR_r2),
                                                               sum_SR_r3 = sum(SR_r3) )

col_sums

```



```{r}

gene_read_CPM = gene_read_counts %>%
    mutate(LR_r1_CPM = LR_r1 / col_sums$sum_LR_r1 * 1e6) %>%
    mutate(LR_r2_CPM = LR_r2 / col_sums$sum_LR_r2 * 1e6) %>%
    mutate(LR_r3_CPM = LR_r3 / col_sums$sum_LR_r3 * 1e6) %>%
    
    mutate(SR_r1_CPM = SR_r1 / col_sums$sum_SR_r1 * 1e6) %>%
    mutate(SR_r2_CPM = SR_r2 / col_sums$sum_SR_r2 * 1e6) %>%
    mutate(SR_r3_CPM = SR_r3 / col_sums$sum_SR_r3 * 1e6) %>%
    
    select(LR_r1_CPM, LR_r2_CPM, LR_r3_CPM,
           SR_r1_CPM, SR_r2_CPM, SR_r3_CPM) 

rownames(gene_read_CPM) = gene_ids

```

```{r}

gene_read_CPM_means = gene_read_CPM %>% rowwise() %>%
    mutate(mean_LR_CPM = mean(LR_r1_CPM, LR_r2_CPM, LR_r3_CPM),
           mean_SR_CPM = mean(SR_r1_CPM, SR_r2_CPM, SR_r3_CPM) )  %>%
    select(mean_LR_CPM, mean_SR_CPM)


gene_read_CPM_means = data.frame(gene_read_CPM_means) 

rownames(gene_read_CPM_means) = gene_ids
```


```{r}

gene_read_CPM_means %>% 
   ggplot(aes(x=mean_LR_CPM, y=mean_SR_CPM)) + geom_point() +
    theme_bw() + 
    scale_x_continuous(trans='log10') +
    scale_y_continuous(trans='log10') +
     stat_smooth(method = "lm", 
              formula = y ~ x, 
              geom = "smooth")


```


```{r}
gene_read_CPM_means_filt = gene_read_CPM_means %>% filter(mean_LR_CPM > 10 & mean_SR_CPM > 10) %>% 
    mutate(SR_div_LR = mean_SR_CPM/mean_LR_CPM) 
```




```{r}

median(gene_read_CPM_means_filt$SR_div_LR)

```

So a median of ~4x as many short reads as long reads per gene


# Examine just the genes involved in the control fusions

```{r}

gene_read_CPM_means_tibble = tibble(gene_read_CPM_means)
gene_read_CPM_means_tibble$gene_id = gene_ids

gene_read_CPM_means_tibble$gene_symbol = sapply(gene_read_CPM_means_tibble$gene_id, function(x) { str_split(x, "\\|")[[1]][1] })

```


```{r}

fusion_gene_cpm_means = NULL

control_fusions = ctatLRF_FI_control_results %>% select(FusionName)  %>% unique() %>% pull(FusionName)

for (control_fusion in control_fusions) {
    
    geneA = str_split(control_fusion, "--")[[1]][1]
    geneB = str_split(control_fusion, "--")[[1]][2]
    
    fusion_gene_cpm_means = bind_rows(fusion_gene_cpm_means,
                                      bind_rows(
                                        gene_read_CPM_means_tibble %>% filter(gene_symbol == geneA),
                                        gene_read_CPM_means_tibble %>% filter(gene_symbol == geneB))  %>% mutate(fusion = control_fusion)  )
}

fusion_gene_cpm_means


```


```{r}

library(ggrepel)

control_fusion_SR_vs_LR_CPM_plot = fusion_gene_cpm_means %>% 
    select(-fusion) %>% unique() %>%
    ggplot(aes(x=mean_LR_CPM, y=mean_SR_CPM)) + geom_point() +
    theme_bw() +
    scale_x_continuous(trans='log10') +
    scale_y_continuous(trans='log10')  +                                      
  stat_smooth(method = "lm", 
              formula = y ~ x, 
              geom = "smooth") +

  geom_abline(slope=1, intercept=0, color='red') +
   geom_text_repel(aes(label=gene_symbol))

control_fusion_SR_vs_LR_CPM_plot

```

```{r}

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

```


```{r}

fusion_gene_cpm_means %>% select(-fusion) %>% unique() %>%
    mutate(SR_div_LR = mean_SR_CPM/mean_LR_CPM) %>%
    summarize(median(SR_div_LR))

```

So even with just the control fusions, we're finding >4x as many short reads as long reads mapping to the corresponding genes.

Note that with short reads we can't determine which non-fusion reads derive from the spiked-in control fusions vs. the normal cell line background.


