---
title: "fitms results smoking ALL"
author: "Keren Xu"
date: "04/28/2022"
output: html_document
---


```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = T)
knitr::opts_chunk$set(warning = FALSE)
knitr::opts_chunk$set(message = FALSE)
knitr::opts_chunk$set(fig.width= 7, fig.height=5) 
```

# Load packages  
```{r}
library(tidyverse)
require(dplyr)
require(data.table)
library(pals)
library(colorblindr)
theme_set(theme_light()) 
setwd("/dir/FitMS.out")
```

## R Markdown
```{r}
# read file path
all_paths <- paste0(list.dirs("/dir/FitMS.out")[-c(1)][seq(3, 152, by = 4)], "/contributions.csv")

# read file content
all_content <-
  all_paths %>%
  lapply(read.csv,
         header = TRUE)
 
# read file name
all_filenames <- all_paths %>%
  as.list()
 
# combine file content list and file name list
all_lists <- mapply(c, all_content, all_filenames, SIMPLIFY = FALSE)
 
# unlist all lists and change column name
all_result <- rbindlist(all_lists, fill = T)

# change column name

all_result$sample_id <- gsub("/signatures/1/contributions.csv", "", all_result$V1) 

all_result$sample_id <- gsub("/dir/FitMS.out/TM_","",all_result$sample_id)

all_result <- all_result %>% select(-V1)

all_result <- all_result %>% filter(Contribution!=0)

all_result %>% count(sample_id)
```

# add unassigned bases 
```{r}
# read file path
all_paths <- paste0(list.dirs("/dir/FitMS.out")[-c(1)][seq(1, 152, by = 4)], "/catalog.csv")

# read file content
all_content <-
  all_paths %>%
  lapply(read.csv,
         header = TRUE)
 
function_sum <- function(df){
  df %>% .[,2] %>% sum()
}
  
total_base <- all_content %>% map(function_sum) %>% unlist()

# read file name
all_filenames <- all_paths 
 
df_total_base <- data.frame(File.Path = all_filenames, total_base = total_base)

df_total_base$sample_id <- gsub("/catalog.csv","", df_total_base$File.Path)

df_total_base$sample_id <- gsub("/dir/FitMS.out/TM_","",df_total_base$sample_id)

df_total_base <- df_total_base %>% select(-c("File.Path"))

df_total_base %>% arrange(-total_base)

```


# stacked bar graph 
```{r}
all_result <- all_result %>% mutate(Contribution = round(Contribution,digits = 0))

all_result <- all_result %>% left_join(all_result %>% group_by(sample_id) %>% summarize(total = sum(Contribution)), by = "sample_id")

all_result <- all_result %>% left_join(df_total_base, by = "sample_id") %>% mutate(unassigned = total_base-total)

unassigned_rows <- all_result %>% select(sample_id, unassigned) %>% distinct() %>% rename("Contribution" = "unassigned") %>% mutate(Signature = "unassigned")

all_result <- rbind(all_result, unassigned_rows, fill = TRUE) %>% select(Signature, Contribution, sample_id)
all_result <- all_result %>% left_join(df_total_base, by = "sample_id") 

# all_result <- all_result %>% mutate(Contribution = ifelse(Contribution<0, 0, Contribution))

all_result %>% count(Signature) %>% pull(Signature) %>% dput()

# Stacked
ggplot(all_result, aes(fill=factor(Signature, levels=c("unassigned", "GEL-Lymphoid_common_SBS127", "GEL-Lymphoid_common_SBS101", "GEL-Lymphoid_common_SBS18", "GEL-Lymphoid_common_SBS17", 
 "GEL-Lymphoid_common_SBS9",  "GEL-Lymphoid_common_SBS5+8", "GEL-Lymphoid_common_SBS5_b", "GEL-Lymphoid_common_SBS5_a", "GEL-Lymphoid_common_SBS1")), y=Contribution, x=reorder(sample_id, -total_base))) + 
    geom_bar(position="stack", stat="identity")+ theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) + scale_fill_manual(values=as.vector(cols25(25))) + guides(fill = guide_legend(title = "Signature"))+labs(x = "sample id")

# compare to the tooth project, it has K but no 
# percent
ggplot(all_result, aes(fill=factor(Signature, levels=c("unassigned", "GEL-Lymphoid_common_SBS127", "GEL-Lymphoid_common_SBS101", "GEL-Lymphoid_common_SBS18", "GEL-Lymphoid_common_SBS17", 
 "GEL-Lymphoid_common_SBS9",  "GEL-Lymphoid_common_SBS5+8", "GEL-Lymphoid_common_SBS5_b", "GEL-Lymphoid_common_SBS5_a", "GEL-Lymphoid_common_SBS1")), y=Contribution, x=reorder(sample_id, -total_base))) + 
    geom_bar(position="fill", stat="identity")+ theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) + scale_fill_manual(values=as.vector(cols25(25))) + guides(fill = guide_legend(title = "Signature"))+labs(x = "sample id")


all_result %>% fwrite("fitms.out.contribution.txt", sep = "\t")
```


# reference contribution
```{r}

# read file path
all_paths <- paste0(list.dirs("/dir/FitMS.out")[-c(1)][seq(3, 152, by = 4)], "/reference_contributions.csv")

# read file content
all_content <-
  all_paths %>%
  lapply(read.csv,
         header = TRUE)
 
# read file name
all_filenames <- all_paths %>%
  as.list()
 
# combine file content list and file name list
all_lists <- mapply(c, all_content, all_filenames, SIMPLIFY = FALSE)
 
# unlist all lists and change column name
all_result <- rbindlist(all_lists, fill = T)

# change column name

all_result$sample_id <- gsub("/signatures/1/reference_contributions.csv", "", all_result$V1) 

all_result$sample_id <- gsub("/dir/FitMS.out/TM_","",all_result$sample_id)

all_result <- all_result %>% select(-V1)

all_result <- all_result %>% filter(Contribution!=0)

all_result %>% count(sample_id)
```


# stacked bar graph 
```{r}
all_result <- all_result %>% mutate(Contribution = round(Contribution,digits = 0))

all_result <- all_result %>% left_join(all_result %>% group_by(sample_id) %>% summarize(total = sum(Contribution)), by = "sample_id")

all_result <- all_result %>% left_join(df_total_base, by = "sample_id") %>% mutate(unassigned = total_base-total)

unassigned_rows <- all_result %>% select(sample_id, unassigned) %>% distinct() %>% rename("Contribution" = "unassigned") %>% mutate(Signature = "unassigned")

all_result <- rbind(all_result, unassigned_rows, fill = TRUE) %>% select(Signature, Contribution, sample_id)
all_result <- all_result %>% left_join(df_total_base, by = "sample_id") 

all_result <- all_result %>% mutate(Contribution = ifelse(Contribution<0, 0, Contribution))


# compare to the tooth project, no RefSig 30, has RefSig N1
# Stacked
ggplot(all_result, aes(fill=factor(Signature, levels=c("unassigned","SBS127","SBS101","SBS18","SBS17","SBS9","SBS8", "SBS5", "SBS1")), y=Contribution, x=reorder(sample_id, -total_base))) + 
    geom_bar(position="stack", stat="identity")+ theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +scale_fill_manual(values=as.vector(cols25(25))) + guides(fill = guide_legend(title = "Signature"))+labs(x = "sample id")


# percent
ggplot(all_result, aes(fill=factor(Signature, levels=c("unassigned","SBS127","SBS101","SBS18","SBS17","SBS9","SBS8", "SBS5", "SBS1")), y=Contribution, x=reorder(sample_id, -total_base))) + 
    geom_bar(position="fill", stat="identity")+ theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) + scale_fill_manual(values=as.vector(cols25(25))) + guides(fill = guide_legend(title = "Signature")) + labs(x = "sample id")


# # check proportions 
# all_result %>% mutate(proportion = Contribution/total_base) -> all_result
# all_result %>% filter(Signature == "RefSig 13") %>% filter(proportion > 0.1) %>% arrange(proportion) %>%  pull(sample_id)
# all_result %>% filter(Signature == "RefSig 2") %>% filter(proportion > 0.1) %>% arrange(proportion) %>%  pull(sample_id)
# all_result %>% filter(Signature == "RefSig 18") %>% filter(proportion > 0.1) %>% arrange(proportion) %>%  pull(sample_id) %>% dput()

all_result %>% fwrite("fitms.out.reference_contribution.txt", sep = "\t")

```


 ˆ