.Rhistory

shiny::runApp()
# Query NYGC data
nygc_junction_data <- query_junction(junc = input$junction_input,
dataset_parquet = "data/nygc/parquets")
runApp()
# Load required libraries
library(shiny)
library(shinydashboard)
library(DT)
library(plotly)
library(tidyverse)
library(arrow)
library(pROC)
library(ggrepel)
library(data.table)
library(glue)
# NYGC plot function
plot_junction_nygc <- function(junc, plotin_table, measure = 'spliced_reads',
gene_name = "") {
vals <- c(`ALS-TDP` = "#E1BE6A", Control = "#40B0A6", `FTLD-TDP` = "#E1BE6A",
`ALS-non-TDP` = "#408A3E", `FTLD-non-TDP` = "#408A3E")
plotin_table <- plotin_table %>%
mutate(tdp_proteinopathy = ifelse(tdp_proteinopathy %in% c("FTLD-TAU","FTLD-FUS"), "FTLD-non-TDP", tdp_proteinopathy))
if(gene_name != ""){
plot_title <- glue::glue("{gene_name} - {junc}")
}else{
plot_title <- glue::glue("{junc}")
}
plotin_table =plotin_table %>%
filter(rna_tissue_source_simplified %in% c("Spinal_Cord_Thoracic", "Spinal_Cord_Cervical", "Spinal_Cord_Lumbar",
"Cortex_Frontal", "Cortex_Motor", "Cortex_Temporal",
"Medulla", "Hippocampus", "Choroid",
"Cortex_Unspecified",
"Spinal_Cord_Unspecified","Cerebellum")) %>%
filter(tdp_proteinopathy != "Unknown")
if(measure == 'spliced_reads') {
plt <- plotin_table %>%
mutate(tdp_proteinopathy = fct_relevel(tdp_proteinopathy, "Control", "ALS-non-TDP", "ALS-TDP", "FTLD-non-TDP", "FTLD-TDP")) %>%
ggplot(aes(x = tdp_proteinopathy, y = junction_count, fill = tdp_proteinopathy)) +
geom_boxplot(outlier.colour = NA) +
geom_jitter(height = 0, alpha = 0.7, pch = 21) +
facet_wrap(~rna_tissue_source_simplified, scales = "free_y") +
scale_fill_manual(values = vals) +
ylab("N spliced reads") +
xlab("") +
ggtitle(plot_title) +
ggpubr::theme_pubr() +
theme(legend.position = 'none') +
theme(text = element_text(size = 10)) +
scale_x_discrete(guide = guide_axis(n.dodge = 2))
} else {
plt <- plotin_table %>%
mutate(tdp_proteinopathy = fct_relevel(tdp_proteinopathy, "Control", "ALS-non-TDP", "ALS-TDP", "FTLD-non-TDP", "FTLD-TDP")) %>%
ggplot(aes(x = tdp_proteinopathy, y = psi, fill = tdp_proteinopathy)) +
geom_boxplot(outlier.colour = NA) +
geom_jitter(height = 0, alpha = 0.7, pch = 21) +
facet_wrap(~rna_tissue_source_simplified, scales = "free_y") +
scale_fill_manual(values = vals) +
ylab("PSI") +
xlab("") +
ggtitle(plot_title) +
ggpubr::theme_pubr() +
theme(legend.position = 'none') +
theme(text = element_text(size = 10)) +
scale_x_discrete(guide = guide_axis(n.dodge = 2)) +
scale_y_continuous(labels = scales::percent_format())
}
return(plt)
}
plot_junction_rimod <- function(junc, table, measure = 'spliced_reads',
gene_name = "") {
vals <- c(`FTD-GRN` = "#E1BE6A", control = "#40B0A6", `FTD-C9` = "#E1BE6A",
`FTD-MAPT` = "#408A3E")
if(gene_name != ""){
plot_title <- glue::glue("{gene_name} - {junc}")
}else{
plot_title <- glue::glue("{junc}")
}
if(measure == 'psi'){
p = table %>%
mutate(DiseaseCode = fct_relevel(DiseaseCode,'control','FTD-MAPT')) %>%
as.data.table() %>%
mutate(psi = junction_count / cluster_count) %>%
ggplot(aes(y = psi,x = DiseaseCode,fill = DiseaseCode)) +
geom_boxplot() +
geom_jitter(height = 0) +
theme_classic() +
ylab("Junction PSI") +
scale_fill_manual(values = vals) +
ylab("PSI") +
xlab("") +
ggtitle(plot_title) +
ggpubr::theme_pubr() +
theme(legend.position = 'none') +
theme(text = element_text(size = 10)) +
scale_y_continuous(labels = scales::percent_format())
}else{
p = table %>%
mutate(DiseaseCode = fct_relevel(DiseaseCode,'control','FTD-MAPT')) %>%
as.data.table() %>%
ggplot(aes(y = junction_count,x = DiseaseCode,fill = DiseaseCode)) +
geom_boxplot() +
geom_jitter(height = 0) +
theme_classic() +
ylab("N spliced reads") +
scale_fill_manual(values = vals) +
xlab("") +
ggtitle(plot_title) +
ggpubr::theme_pubr() +
theme(legend.position = 'none') +
theme(text = element_text(size = 10))
}
return(p)
}
# query the parquets files
query_junction = function(junc, dataset_parquet){
parts <- unlist(strsplit(junc, "[:\\-]"))
queried_chrom = parts[[1]]
queried_start = as.numeric(parts[[2]])
queried_end = as.numeric(parts[[3]])
dataset_parquet_query = arrow::open_dataset(dataset_parquet)
junction_data <- dataset_parquet_query %>%
filter(chrom == queried_chrom,
start == queried_start,
end == queried_end) %>%
collect()
return(junction_data)
}
# Query NYGC data
nygc_junction_data <- query_junction(junc = input$junction_input,
dataset_parquet = "data/nygc/parquets")
input = list()
input$junction_input
input$junction_input = chr8:79611214-79636802
input$junction_input = "chr8:79611214-79636802"
# Query NYGC data
nygc_junction_data <- query_junction(junc = input$junction_input,
dataset_parquet = "data/nygc/parquets")
nygc_junction_data
runApp()
nygc_junction_data = nygc_junction_data %>%
distinct(chrom,sample,cluster_count,start,end,cluster,strand,junction_count,psi)
nygc_junction_data
# Query RIMOD data
rimod_junction_data <- query_junction(junc = input$junction_input,
dataset_parquet = "data/rimod/parquets")
rimod_junction_data
runApp()
runApp()
runApp()
runApp()
runApp()
possible_events_dt <- fread("data/rimod/rimod_meta.csv")
possible_events_dt
runApp()
runApp()
# If we've got a gene name for that event in this table
junction_name = possible_events_dt %>%
filter(junctions_coords == input$junction_input) %>%
distinct(gene_name) %>%
pull(gene_name)
junction_name
nchar(junction_name)
if(nchar(junction_name) == 0){
junction_name = ""
}
if(nchar(junction_name) == 0){
junction_name = ""
}
nchar(junction_name) == 0
nchar(junction_name) == 0
junction_name
runApp()
junction_name
nchar(junction_name)
nchar(junction_name) == 0
possible_events_dt %>%
filter(junctions_coords == input$junction_input) %>%
distinct(gene_name)
possible_events_dt %>%
filter(junctions_coords == input$junction_input) %>%
distinct(gene_name) %>% nrow()
row_num = possible_events_dt %>%
filter(junctions_coords == input$junction_input) %>%
distinct(gene_name) %>% nrow()
possible_events_dt %>% row_num
row_num
runApp()
shiny::runApp()
shiny::runApp()
runApp()
runApp()
runApp()
runApp()
runApp()
runApp()
runApp()
runApp()
runApp()
runApp()
runApp()
runApp()
runApp()
runApp()
shiny::runApp()
library(tximeta)
library(tidyverse)
library(fishpond)
library(tximport)
library(ggpubr)
library(data.table)
library(SummarizedExperiment)
mapped_size = fread("~/cluster/vyplab/first_weeks/recalled_nanopore/aligned_library_size.csv")
skey = fread("/Users/annaleigh/Documents/GitHub/m6a_analyses/samplekey.csv",fill = TRUE)
mapped_size = mapped_size %>%
mutate(genome = ifelse(grepl("New_UNC13A_MiniGene_2xH",sample_id),'mini_gene',"genome")) %>%
separate(sample_id,into = 'sample_id') %>%
left_join(skey %>% select(sample_id,sample_name)) %>%
separate(sample_name,into = c("name"))  %>%
mutate(condition = ifelse(grepl("Knockdown",name),"TDP43_KD","Control")) %>%
mutate(library_size_mr = library_size / 10^6)
m = list.files('~/cluster/vyplab/first_weeks/recalled_nanopore/oarfish/',pattern = '.quant',full.names = TRUE)
gsub(".quant","",gsub(".*//","",m))
names(m) = mapped_size$name[match(gsub(".quant","",gsub(".*//","",m)), mapped_size$sample_id)]
tx2gene = fread("tx2gene_oarfish.csv")
tx2gene = tx2gene %>%
mutate(gene_id = gsub("\\..*","",geneCol)) %>%
select(-geneCol)
k_skip = tximeta(coldata = data.frame(files = m,
names = names(m)),
type = 'oarfish',
skipMeta = TRUE)