---
title: "M132TS_analysis"
author: "bhaas"
date: '2023-09-28'
output: github_document
---

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


```{r}
MIN_CELLS = 5
```


# melanoma tumor sample M132TS 


```{r}

# parse fusion read support info including cell barcodes and umis

Tum_data = read.table("data/M132TS.melanoma_sc.filtered_cells_and_dedup_umis.tsv.gz", header=T, sep="\t", stringsAsFactors = F)%>%
    mutate(celltype_final = ifelse(leiden %in% c(4,7), 'tumor', 'normal'))

head(Tum_data)


```


```{r}
# since starF and FI were run at max sensitivity, lets restrict fusions to those identified by ctat-LRF

ctat_LRF_fusion_genes = Tum_data %>% filter(method == 'ctat-LR-fusion') %>% select(FusionName) %>% unique() %>% pull(FusionName)

Tum_data = Tum_data %>% filter(FusionName %in% ctat_LRF_fusion_genes)
```


```{r}

fusion_annots = read.table("data/M132TS.fusion_annots.gz", sep="\t", header=T, stringsAsFactors = F)

```


```{r}

Tum_umap_data = read.table("data/M132TS.bc_to_umap_n_leiden.tsv", header=T, sep="\t", stringsAsFactors = F) 

# number of cells
num_tumor_cells = nrow(Tum_umap_data) 
message("number tumor cells: ", num_tumor_cells)
# 6932 total cells

Tum_umap_data %>% group_by(leiden) %>% tally(name='count_cell_cluster') %>% mutate(frac_tot_cells = count_cell_cluster/num_tumor_cells)
# cluster ids 4,7 = tumor
# 584 + 117 = 701 tumor cells


```
```{r}

Tum_cell_counts = Tum_data %>% select(FusionName, cell_barcode) %>% unique() %>% 
    group_by(FusionName) %>%
    tally(name='tot_cells_w_fusion') %>%
    mutate(frac_tot_cells=tot_cells_w_fusion/num_tumor_cells)  %>%
    arrange(desc(frac_tot_cells))

Tum_cell_counts = left_join(Tum_cell_counts, fusion_annots)

Tum_cell_counts %>% filter(tot_cells_w_fusion >= MIN_CELLS)
```

```{r}
fusions_min_cell_counts = Tum_cell_counts %>% filter(tot_cells_w_fusion >= MIN_CELLS) 

fusions_min_cell_counts

```





```{r}

# examine distribution of fusion calls according to cell types

Tum_fusion_frac_cell_types = Tum_data %>% select(FusionName, cell_barcode, celltype_final) %>% unique() %>%
    group_by(FusionName, celltype_final) %>% tally(name='tot_cells_w_fusion') %>% 
    mutate(frac_fusion_cells=prop.table(tot_cells_w_fusion)) %>%
    arrange(desc(tot_cells_w_fusion))

Tum_fusion_frac_cell_types %>% head()


```


```{r}

Tum_data %>% select(method) %>% unique()

```




```{r}

starF_fusions = Tum_data %>% filter(method=="STAR-Fusion")

FI_fusions = Tum_data %>% filter(method=="FusionInspector")

ctat_LRF_fusions = Tum_data %>% filter(method == "ctat-LR-fusion")


```






```{r}

# parse cell counts by detection method

Tum_cell_counts_by_method = read.table("data/M132TS.melanoma_sc.filtered_cells_and_dedup_umis.cell_counts_by_method.tsv.gz",
                                          header=T, sep="\t", stringsAsFactors = F)


Tum_cell_counts_by_method  = Tum_cell_counts_by_method %>% filter(FusionName %in% ctat_LRF_fusion_genes)

Tum_cell_counts_by_method %>% head()

```

```{r}

# reorganize to compare findings of cells according to fusion and method

Tum_cell_counts_by_method_spread = Tum_cell_counts_by_method %>% select(FusionName, LeftBreakpoint, RightBreakpoint, method, leiden, cell_counts) %>%
    spread(key=method, value=cell_counts) %>% 
    arrange(desc(`ctat-LR-fusion`)) 

Tum_cell_counts_by_method_spread %>% head()

```



# plot counts of cells for these fusions:

```{r}

right_join(Tum_cell_counts_by_method, 
          Tum_cell_counts_by_method %>% 
                         filter(cell_counts >= MIN_CELLS)  %>% 
                         select(FusionName, LeftBreakpoint, RightBreakpoint) 
          ) %>%
              rowwise() %>% mutate(fusion=paste(FusionName, LeftBreakpoint, RightBreakpoint, collapse=":")) %>%
              ggplot(aes(x=fusion, y=cell_counts, fill=method)) + geom_bar(stat='identity', position='dodge') +
              theme(axis.text.x = element_text(angle = 90, hjust = 1))


```



# Examine cell type representation by fusions


```{r}

fusion_frac_cell_types = Tum_data %>% select(FusionName, cell_barcode, leiden) %>% unique() %>%
    group_by(FusionName, leiden) %>% tally(name='tot_cells_w_fusion') %>% 
    mutate(frac_fusion_cells=prop.table(tot_cells_w_fusion)) %>%
    arrange(desc(tot_cells_w_fusion))

fusion_frac_cell_types %>% head()

```




```{r}

fusions_of_interest = Tum_fusion_frac_cell_types %>% 
    filter(FusionName %in% fusions_min_cell_counts$FusionName) %>%
    arrange(desc(tot_cells_w_fusion)) %>%
    filter(celltype_final == "tumor" & frac_fusion_cells >= 0.8)

fusions_of_interest

```


```{r}

fusions_of_interest = left_join(fusions_of_interest,
          fusion_annots)

fusions_of_interest
```


NUTM2A-AS1--RP11-203L2.4 is the only relevant inter-chromosomal.

The others have few cells and are likely cis-spliced fusions.




```{r}

baseplot = Tum_umap_data %>% ggplot(aes(x=umap_1, y=umap_2)) + theme_bw() + geom_point()

baseplot + geom_point(aes(color=as.factor(leiden)))

```


```{r}

Tum_umap_data %>% ggplot(aes(x=umap_1, y=umap_2)) + theme_bw() + geom_point(aes(color=CST3))

```

```{r}

TUMOR_CLUSTERS = c(4,7)

Tum_umap_data = Tum_umap_data %>% mutate(tumor = leiden %in% TUMOR_CLUSTERS)

Tum_umap_data %>%
    ggplot(aes(x=umap_1, y=umap_2)) + theme_bw() + geom_point(aes(color=tumor))


```



```{r}



Tum_umap_data %>% select(index, tumor) %>% unique() %>% nrow()

# 6932 total cells
    
 Tum_umap_data %>% select(index, tumor) %>% unique() %>% group_by(tumor) %>% tally()
 
# tumor #cells
# FALSE	6231			
# TRUE	701	

NUM_TOTAL_CELLS = 6932
NUM_TUMOR_CELLS = 701
NUM_NORMAL_CELLS = 6231

# so 10% are tumor cells
```




```{r}

# label tumor cells
Tum_data = Tum_data %>% 
    mutate(tumor = leiden %in% TUMOR_CLUSTERS)

fusion_frac_tumor = Tum_data %>%
    select(FusionName, cell_barcode, tumor) %>% unique() %>%
    group_by(FusionName, tumor) %>% tally(name='tot_cells_w_fusion') %>%
    spread(key=tumor, value=tot_cells_w_fusion, fill = 0) %>%
    rename(normal_cell_count=`FALSE`) %>% rename(tumor_cell_count=`TRUE`) %>%
    mutate(frac_normal_cells = normal_cell_count / NUM_NORMAL_CELLS,
           frac_tumor_cells = tumor_cell_count / NUM_TUMOR_CELLS) %>%
    arrange(desc(frac_tumor_cells))



fusions_of_interest = fusion_frac_tumor %>% filter(frac_tumor_cells > 0.01 | frac_normal_cells > 0.01)

fusions_of_interest

# Only one fusion was found in more than 1% of tumor or normal cells: NUTM2A-AS1::RP11-203L2.4 found in 265 tumor cells (38%) and only 3 normal cells (0.05%) through a combination of long and short read fusion transcript analyses (Figure 7a,b);

```




```{r}

x = 0

plots = list()

for (fusion in  fusions_of_interest$FusionName) {
    
    p = baseplot + geom_point(data=Tum_data %>% filter(FusionName == fusion) %>% select(umap_1, umap_2, method) %>% unique(), 
                              aes(color=method), alpha=0.5, size=rel(2)) + 
        ggtitle(paste("M132TS, Fusion: ", fusion) )
    
    plot(p)   
    
    x = x+1
    plots[[x]] = p
}



```




```{r}


pdf("M132TS.fusions_of_interest.pdf")
for (p in plots) {
    plot(p)
}

dev.off()

```


```{r}

fusions_of_interest = left_join(fusions_of_interest,
          fusion_annots)

fusions_of_interest

```

```{r}

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


```



```{r}

Tum_cell_counts_by_method %>% filter(FusionName %in% fusions_of_interest$FusionName) %>% 
    select(FusionName, LeftBreakpoint, RightBreakpoint, method, leiden, cell_counts) %>%
    spread(key=method, value=cell_counts) %>% 
    arrange(desc(`ctat-LR-fusion`))

```



```{r}

M132TS_fusions_of_interest_counts = Tum_data %>% filter(FusionName %in% fusions_of_interest$FusionName) %>% 
        filter(tumor) %>%
        select(FusionName, method, cell_barcode) %>% unique() %>%
        group_by(FusionName, method) %>% tally(name='cell_counts')

M132TS_fusions_of_interest_counts

M132TS_fusions_of_interest_counts %>%
              ggplot(aes(x=FusionName, y=cell_counts, fill=method)) + geom_bar(stat='identity', position='dodge') +
              theme(axis.text.x = element_text(angle = 90, hjust = 1)) +
    ggtitle("M132TS Fusions of Interest: Cell Counts")


```


# focus on NUTM2A-AS1--RP11-203L2.4


```{r}

Tum_data %>% filter(FusionName == "NUTM2A-AS1--RP11-203L2.4" ) %>% 
        filter(tumor) %>%
        select(FusionName, method, cell_barcode) %>% unique() %>%
        group_by(FusionName, method) %>% tally(name='cell_counts') %>%
              ggplot(aes(x=FusionName, y=cell_counts, fill=method)) + geom_bar(stat='identity', position='dodge') +
              theme(axis.text.x = element_text(angle = 90, hjust = 1)) +
    ggtitle("M132TS Fusions of Interest: Cell Counts")


```

```{r}
Tum_data %>% filter(FusionName == "NUTM2A-AS1--RP11-203L2.4" ) %>% 
        filter(tumor) %>%
        select(FusionName, method, cell_barcode) %>% unique() %>%
        group_by(FusionName, method) %>% tally(name='cell_counts') %>%
    mutate(frac_NUTM2A_AS1_fusion_positive = cell_counts / 265)

```






```{r}

 p = baseplot + geom_point(data=Tum_data %>% filter(FusionName == "NUTM2A-AS1--RP11-203L2.4") %>% select(umap_1, umap_2, method) %>% unique(), 
                              aes(color=method), alpha=0.5, size=rel(2)) + 
        ggtitle(paste("M132TS, Fusion: ", fusion) )
    
    plot(p)   



```

```{r}

ggsave(p, file="M132TS_NUTM2A-AS1_fusion_umap.svg", width=7, height=4)


```



```{r}

venn_dist = Tum_data %>% filter(FusionName == "NUTM2A-AS1--RP11-203L2.4") %>%
    filter(tumor) %>%
    select(cell_barcode, method) %>% unique() %>%
    group_by(cell_barcode) %>% arrange(method) %>%
    mutate(methods = paste(method, collapse=",")) %>%
    ungroup() %>%
    select(cell_barcode, methods) %>% unique() %>%
    group_by(methods) %>% tally(name='num_cells')

venn_dist = venn_dist %>% mutate(frac_pos_cells = num_cells/ 265)

venn_dist 

# ~60% of cells identified by ctat-LRF (long reads) alone.

```
60% identified by ctat-LRF only.


```{r}

# what fraction of fusion positive cells were by short reads only?
venn_dist %>% filter(! grepl("ctat-LR-fusion", methods))

venn_dist %>% filter(! grepl("ctat-LR-fusion", methods)) %>% summarize(sum(frac_pos_cells))

# ~ 20% uniquely by short reads only

# and so the remaining ~20% invvolved a combo of short and long reads

```
~20% by short reads alone


```{r}

# double check number of tumor cells with NUTM2A-AS1 fusion = 265

venn_dist %>% summarize(sum(num_cells))


```


# Other fusions of interest in normal cells?

```{r}
# examine how these fusions are distributed among cell clusters



Tum_data %>% 
    filter(FusionName %in% fusions_min_cell_counts) %>%
    select(FusionName, cell_barcode, celltype_final) %>% 
    unique() %>%
    group_by(FusionName, celltype_final) %>% tally(name='fusion_cell_counts_per_cluster') %>% 
    mutate(frac_fusion_cells=prop.table(fusion_cell_counts_per_cluster)) %>%
    arrange(FusionName, desc(fusion_cell_counts_per_cluster)) 


```


```{r}
# just top cell type

# examine how these fusions are distributed among cell clusters

left_join(Tum_data %>% filter(FusionName %in% fusions_min_cell_counts$FusionName) %>% select(FusionName, cell_barcode, leiden) %>% 
    unique() %>%
    group_by(FusionName, leiden) %>% tally(name='fusion_cell_counts_per_cluster') %>% ungroup() %>%
    mutate(frac_fusion_cells=prop.table(fusion_cell_counts_per_cluster)) %>%
    arrange(FusionName, desc(fusion_cell_counts_per_cluster))  %>%
    mutate(tumor_or_normal = (leiden %in% c(4,7))) %>%
    group_by(FusionName) %>% arrange(desc(fusion_cell_counts_per_cluster)) %>% filter(row_number()==1) %>% ungroup() %>%
        arrange(desc(fusion_cell_counts_per_cluster)),
    
    fusion_annots, by='FusionName')


```

# exploring individual fusions 


```{r}

report_on_fusion = function(fusion_name) {
    
    print(Tum_cell_counts_by_method_spread %>% filter(FusionName == fusion_name))
    
    print(Tum_cell_counts %>% filter(FusionName == fusion_name) %>% mutate(type="Tum"))
    
    print(fusion_frac_cell_types %>% filter(FusionName == fusion_name) %>% mutate(type="Tum"))
    
    print(Tum_fusion_frac_cell_types %>% filter(FusionName == fusion_name) %>% mutate(type="Tum"))
    
}


```

```{r}

report_on_fusion("RP11-444D3.1--SOX5")

# found in 7 malignant and 3 normal cells 

```


10 cells with RP11-444D3.1--SOX5 fusion, 70% are in tumor, 30% are in normal clusters.

## RP11-208G20.2--PSPHP1

```{r}

report_on_fusion("RP11-208G20.2--PSPHP1")

```

none here.




```{r}

report_on_fusion("RP1-34H18.1--NAV3")

```






