---
title: "TMM_Muscle"
output: html_document
date: '2023-12-04'
---
```{r}
library()
```

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```
```{r message=FALSE, warning=FALSE, paged.print=FALSE}
library(arrow,lib.loc = '~/userRLib/')
library(cowplot)
library(dplyr)
library(tidyr)
library(PCAtools)
library(boot)
library(edgeR)
library(tibble)
library(stringr)
library(ggplot2)
library(optparse)
library(GenomicFeatures)
```
```{r}
# methods
object <- '/oak/stanford/groups/smontgom/tannerj/UDN_LRS_watershed/Watershed-SV/new-protein-lincRNA-100k-rare-position-aware-2.0/eval_watarshed_ABC_regional_tissue_tmm_corrected__SV_Watershed_approximate_100000_number_of_dimensions_3_pthresh_0.0027_pval_frac_0.098_shared_trainP_2.56__evaluation_object.rds'
data <- '/oak/stanford/groups/smontgom/tannerj/UDN_LRS_watershed/Watershed-SV/new-protein-lincRNA-100k-rare-position-aware-2.0/eval_watarshed_ABC_regional.no_vep.z3.tsv'

# compute aucprc

remove_na <- function(x) {
	return(x[!is.na(x)])
}

compute_auc_across_dimensions <- function(number_of_dimensions, watershed_posteriors, gam_posteriors, real_valued_outliers_test1, binary_outliers_test2) {
	# Initialize objects to store auc information
	auc_object_across_dimensions <- list()
	pos_list <- c()
	neg_list <- c()
	for (dimension in 1:number_of_dimensions) {
  		# Pseudo gold standard
  		test_outlier_status <- binary_outliers_test2[,dimension]
  		# Watershed evaluation curves
  		watershed_roc_obj <- roc.curve(scores.class0 = remove_na(watershed_posteriors[,dimension][test_outlier_status==1 & !is.na(real_valued_outliers_test1[,dimension])]), scores.class1 = remove_na(watershed_posteriors[,dimension][test_outlier_status==0 & !is.na(real_valued_outliers_test1[,dimension])]), curve = T)
  		watershed_pr_obj <- pr.curve(scores.class0 = remove_na(watershed_posteriors[,dimension][test_outlier_status==1 & !is.na(real_valued_outliers_test1[,dimension])]), scores.class1 = remove_na(watershed_posteriors[,dimension][test_outlier_status==0 & !is.na(real_valued_outliers_test1[,dimension])]), curve = T)
  	
  		# GAM evaluation curves
  		gam_roc_obj <- roc.curve(scores.class0 = remove_na(gam_posteriors[,dimension][test_outlier_status==1 & !is.na(real_valued_outliers_test1[,dimension])]), scores.class1 = remove_na(gam_posteriors[,dimension][test_outlier_status==0 & !is.na(real_valued_outliers_test1[,dimension])]), curve = T)
   		gam_pr_obj <- pr.curve(scores.class0 = remove_na(gam_posteriors[,dimension][test_outlier_status==1 & !is.na(real_valued_outliers_test1[,dimension])]), scores.class1 = remove_na(gam_posteriors[,dimension][test_outlier_status==0 & !is.na(real_valued_outliers_test1[,dimension])]), curve = T)

		evaROC <-	
		 list(watershed_sens=watershed_roc_obj$curve[,2],
              watershed_spec=1-watershed_roc_obj$curve[,1],
         	  watershed_auc=watershed_roc_obj$auc,
         	  watershed_pr_auc=watershed_pr_obj$auc.integral,
         	  watershed_recall=watershed_pr_obj$curve[,1],
         	  watershed_precision=watershed_pr_obj$curve[,2],
         	  GAM_sens=gam_roc_obj$curve[,2],
              GAM_spec=1-gam_roc_obj$curve[,1],
              GAM_auc=gam_roc_obj$auc,
         	  GAM_pr_auc=gam_pr_obj$auc.integral,
         	  GAM_recall=gam_pr_obj$curve[,1],
         	  GAM_precision=gam_pr_obj$curve[,2])

		auc_object_across_dimensions[[dimension]] <- list(evaROC=evaROC)
	}
	return(auc_object_across_dimensions)
}
# wrapper function to fit compute_auc_pr_across_dimensions into boot workflow. 
stat_compute_delta_auprc <- function(eval_df,sampled_indices,tissues=1) {
  watershed <- matrix(eval_df$watershed_n2_posterior[sampled_indices,1])
  gam <- matrix(eval_df$gam_n2_posterior[sampled_indices,1])
  rvot1 <- matrix(eval_df$rvot1[sampled_indices,1])
  bot2 <- matrix(eval_df$bot2[sampled_indices,1])
  auc_object_across_dimensions <- compute_auc_across_dimensions(tissues,
                                                              matrix(watershed),
                                                              matrix(gam), 
                                                              matrix(rvot1), 
                                                              matrix(bot2))
  return (auc_object_across_dimensions[[1]]$evaROC$watershed_pr_auc - auc_object_across_dimensions[[1]]$evaROC$GAM_pr_auc)

}

# calculate auprc
re_calculate_precision_recall <- function(watershed_object,eval_data,tissue_count,fraction_outlier,pvalue_threshold,compute_ci = FALSE) {
  watershed_ABC_100k_obj <-readRDS(watershed_object)
  data_input <- load_watershed_data(eval_data, tissue_count, fraction_outlier, pvalue_threshold)

  # Parse data_input for evaluation-related relevant fields
  feat_all <- data_input$feat
  discrete_outliers_all <- data_input$outliers_discrete
  binary_outliers_all <- data_input$outliers_binary
  fraction_binary_outliers_all <- data_input$fraction_outliers_binary
  N2_pairs <- data_input$N2_pairs
  
  # Extraction of test data for Watershed
  if (tissue_count!=1) {
    feat_test <- rbind(feat_all[(!is.na(N2_pairs)),][seq(from=1,to=sum((!is.na(N2_pairs))),by=2),], 
                       feat_all[(!is.na(N2_pairs)),][seq(from=2,to=sum((!is.na(N2_pairs))),by=2),])
    
    discrete_outliers_test1 <- rbind(discrete_outliers_all[!is.na(N2_pairs),][seq(from=1,to=sum(!is.na(N2_pairs)),by=2),],
                                     discrete_outliers_all[!is.na(N2_pairs),][seq(from=2,to=sum(!is.na(N2_pairs)),by=2),])
    discrete_outliers_test2 <- rbind(discrete_outliers_all[!is.na(N2_pairs),][seq(from=2,to=sum(!is.na(N2_pairs)),by=2),],
                                     discrete_outliers_all[!is.na(N2_pairs),][seq(from=1,to=sum(!is.na(N2_pairs)),by=2),])
    binary_outliers_test1 <- rbind(binary_outliers_all[!is.na(N2_pairs),][seq(from=1,to=sum(!is.na(N2_pairs)),by=2),],
                                   binary_outliers_all[!is.na(N2_pairs),][seq(from=2,to=sum(!is.na(N2_pairs)),by=2),])
    binary_outliers_test2 <- rbind(fraction_binary_outliers_all[!is.na(N2_pairs),][seq(from=2,to=sum(!is.na(N2_pairs)),by=2),],
                                   fraction_binary_outliers_all[!is.na(N2_pairs),][seq(from=1,to=sum(!is.na(N2_pairs)),by=2),])
    binary_outliers_test2_alt <- rbind(binary_outliers_all[!is.na(N2_pairs),][seq(from=2,to=sum(!is.na(N2_pairs)),by=2),],
                                   binary_outliers_all[!is.na(N2_pairs),][seq(from=1,to=sum(!is.na(N2_pairs)),by=2),])
    real_valued_outliers_test1 <- -log10(abs(rbind(data_input$outlier_pvalues[!is.na(N2_pairs),][seq(from=1,to=sum(!is.na(N2_pairs)),by=2),],
                                                   data_input$outlier_pvalues[!is.na(N2_pairs),][seq(from=2,to=sum(!is.na(N2_pairs)),by=2),])) + 1e-7)}
  else {
    feat_test <- rbind(feat_all[!is.na(N2_pairs),][seq(from=1,to=sum(!is.na(N2_pairs)),by=2),], feat_all[!is.na(N2_pairs),][seq(from=2,to=sum(!is.na(N2_pairs)),by=2),])
    discrete_outliers_test1 <- as.matrix(c(discrete_outliers_all[!is.na(N2_pairs)][seq(from=1,to=sum(!is.na(N2_pairs)),by=2)],
                                           discrete_outliers_all[!is.na(N2_pairs)][seq(from=2,to=sum(!is.na(N2_pairs)),by=2)]))
    
    discrete_outliers_test2 <- as.matrix(c(discrete_outliers_all[!is.na(N2_pairs)][seq(from=2,to=sum(!is.na(N2_pairs)),by=2)],
                                           discrete_outliers_all[!is.na(N2_pairs)][seq(from=1,to=sum(!is.na(N2_pairs)),by=2)]))
    
    binary_outliers_test1 <- as.matrix(c(binary_outliers_all[!is.na(N2_pairs)][seq(from=1,to=sum(!is.na(N2_pairs)),by=2)],
                                         binary_outliers_all[!is.na(N2_pairs)][seq(from=2,to=sum(!is.na(N2_pairs)),by=2)]))
    
    binary_outliers_test2 <- as.matrix(c(fraction_binary_outliers_all[!is.na(N2_pairs)][seq(from=2,to=sum(!is.na(N2_pairs)),by=2)],
                                         fraction_binary_outliers_all[!is.na(N2_pairs)][seq(from=1,to=sum(!is.na(N2_pairs)),by=2)]))
    # Absolute pvalues from test prediction data set (to be used for RNA-only analysis)
    real_valued_outliers_test1 <- -log10(abs(as.matrix(c(data_input$outlier_pvalues[!is.na(N2_pairs)][seq(from=1,to=sum(!is.na(N2_pairs)),by=2)],
                                                         data_input$outlier_pvalues[!is.na(N2_pairs)][seq(from=2,to=sum(!is.na(N2_pairs)),by=2)]))) + 1e-7)
    
    binary_outliers_test2_alt <- as.matrix(c(binary_outliers_all[!is.na(N2_pairs)][seq(from=2,to=sum(!is.na(N2_pairs)),by=2)],
                                             binary_outliers_all[!is.na(N2_pairs)][seq(from=1,to=sum(!is.na(N2_pairs)),by=2)]))

  }
  # recalculate AUC
  if ( compute_ci ==TRUE ) {
    eval_data <- data_frame(watershed_n2_posterior = watershed_ABC_100k_obj$n2_posteriors$watershed_n2_posterior,
                            gam_n2_posterior = watershed_ABC_100k_obj$n2_posteriors$gam_n2_posterior,
                            rvot1 = real_valued_outliers_test1,
                            bot2 = binary_outliers_test2_alt)
    bootstrap_delta_auprc <- boot(eval_data,statistic = stat_compute_delta_auprc, R = 5000)
    bootstrap_ci_delta_auprc <- boot.ci(bootstrap_delta_auprc)
    return (list(boot_out = bootstrap_delta_auprc, boot_ci = bootstrap_ci_delta_auprc))
  } else {
  watershed <- watershed_ABC_100k_obj$n2_posteriors$watershed_n2_posterior
  gam <- watershed_ABC_100k_obj$n2_posteriors$gam_n2_posterior
  rvot1 <- real_valued_outliers_test1
  bot2 <- binary_outliers_test2_alt
  auc_object_across_dimensions <- compute_auc_across_dimensions(1,
                                                                watershed,
                                                                gam, 
                                                                rvot1, 
                                                                bot2)
  return (auc_object_across_dimensions)
  }
}


re_calculate_precision_recall_noncoding <- function(watershed_object,eval_data,tissue_count,fraction_outlier,pvalue_threshold,compute_ci = FALSE) {
  watershed_ABC_100k_obj <-readRDS(watershed_object)
  data_input <- load_watershed_data(eval_data, tissue_count, fraction_outlier, pvalue_threshold)

  # Parse data_input for evaluation-related relevant fields
  feat_all <- data_input$feat
  discrete_outliers_all <- data_input$outliers_discrete
  binary_outliers_all <- data_input$outliers_binary
  fraction_binary_outliers_all <- data_input$fraction_outliers_binary
  N2_pairs <- data_input$N2_pairs
  
  # Extraction of test data for Watershed
  if (tissue_count!=1) {
    feat_test <- feat_all[(!is.na(N2_pairs)),]
    is_noncoding_1 <- c((feat_test[,"coding_fraction_affected" ]==0)[seq(from=1,to=sum((!is.na(N2_pairs))),by=2)],
                  (feat_test[,"coding_fraction_affected" ]==0)[seq(from=2,to=sum((!is.na(N2_pairs))),by=2)])
    is_noncoding_2 <- c((feat_test[,"coding_fraction_affected" ]==0)[seq(from=2,to=sum((!is.na(N2_pairs))),by=2)],
                  (feat_test[,"coding_fraction_affected" ]==0)[seq(from=1,to=sum((!is.na(N2_pairs))),by=2)])

    discrete_outliers_test1 <- rbind(discrete_outliers_all[!is.na(N2_pairs),][seq(from=1,to=sum(!is.na(N2_pairs)),by=2),],
                                     discrete_outliers_all[!is.na(N2_pairs),][seq(from=2,to=sum(!is.na(N2_pairs)),by=2),])
    discrete_outliers_test2 <- rbind(discrete_outliers_all[!is.na(N2_pairs),][seq(from=2,to=sum(!is.na(N2_pairs)),by=2),],
                                     discrete_outliers_all[!is.na(N2_pairs),][seq(from=1,to=sum(!is.na(N2_pairs)),by=2),])
    binary_outliers_test1 <- rbind(binary_outliers_all[!is.na(N2_pairs),][seq(from=1,to=sum(!is.na(N2_pairs)),by=2),],
                                   binary_outliers_all[!is.na(N2_pairs),][seq(from=2,to=sum(!is.na(N2_pairs)),by=2),])
    binary_outliers_test2 <- rbind(fraction_binary_outliers_all[!is.na(N2_pairs),][seq(from=2,to=sum(!is.na(N2_pairs)),by=2),],
                                   fraction_binary_outliers_all[!is.na(N2_pairs),][seq(from=1,to=sum(!is.na(N2_pairs)),by=2),])
    binary_outliers_test2_alt <- rbind(binary_outliers_all[!is.na(N2_pairs),][seq(from=2,to=sum(!is.na(N2_pairs)),by=2),],
                                   binary_outliers_all[!is.na(N2_pairs),][seq(from=1,to=sum(!is.na(N2_pairs)),by=2),])
    real_valued_outliers_test1 <- -log10(abs(rbind(data_input$outlier_pvalues[!is.na(N2_pairs),][seq(from=1,to=sum(!is.na(N2_pairs)),by=2),],
                                                   data_input$outlier_pvalues[!is.na(N2_pairs),][seq(from=2,to=sum(!is.na(N2_pairs)),by=2),])) + 1e-7)}
  else {
    feat_test <- feat_all[!is.na(N2_pairs),]
    is_noncoding_1 <- c((feat_test[,"coding_fraction_affected" ]==0)[seq(from=1,to=sum((!is.na(N2_pairs))),by=2)],
                  (feat_test[,"coding_fraction_affected" ]==0)[seq(from=2,to=sum((!is.na(N2_pairs))),by=2)])
    is_noncoding_2 <- c((feat_test[,"coding_fraction_affected" ]==0)[seq(from=2,to=sum((!is.na(N2_pairs))),by=2)],
                  (feat_test[,"coding_fraction_affected" ]==0)[seq(from=1,to=sum((!is.na(N2_pairs))),by=2)])

    discrete_outliers_test1 <- as.matrix(c(discrete_outliers_all[!is.na(N2_pairs)][seq(from=1,to=sum(!is.na(N2_pairs)),by=2)],
                                           discrete_outliers_all[!is.na(N2_pairs)][seq(from=2,to=sum(!is.na(N2_pairs)),by=2)]))
    
    discrete_outliers_test2 <- as.matrix(c(discrete_outliers_all[!is.na(N2_pairs)][seq(from=2,to=sum(!is.na(N2_pairs)),by=2)],
                                           discrete_outliers_all[!is.na(N2_pairs)][seq(from=1,to=sum(!is.na(N2_pairs)),by=2)]))
    
    binary_outliers_test1 <- as.matrix(c(binary_outliers_all[!is.na(N2_pairs)][seq(from=1,to=sum(!is.na(N2_pairs)),by=2)],
                                         binary_outliers_all[!is.na(N2_pairs)][seq(from=2,to=sum(!is.na(N2_pairs)),by=2)]))
    
    binary_outliers_test2 <- as.matrix(c(fraction_binary_outliers_all[!is.na(N2_pairs)][seq(from=2,to=sum(!is.na(N2_pairs)),by=2)],
                                         fraction_binary_outliers_all[!is.na(N2_pairs)][seq(from=1,to=sum(!is.na(N2_pairs)),by=2)]))
    # Absolute pvalues from test prediction data set (to be used for RNA-only analysis)
    real_valued_outliers_test1 <- -log10(abs(as.matrix(c(data_input$outlier_pvalues[!is.na(N2_pairs)][seq(from=1,to=sum(!is.na(N2_pairs)),by=2)],
                                                         data_input$outlier_pvalues[!is.na(N2_pairs)][seq(from=2,to=sum(!is.na(N2_pairs)),by=2)]))) + 1e-7)
    
    binary_outliers_test2_alt <- as.matrix(c(binary_outliers_all[!is.na(N2_pairs)][seq(from=2,to=sum(!is.na(N2_pairs)),by=2)],
                                             binary_outliers_all[!is.na(N2_pairs)][seq(from=1,to=sum(!is.na(N2_pairs)),by=2)]))

  }
  if ( compute_ci ==TRUE ) {
    # recalculate AUC
    eval_data <- data_frame(watershed_n2_posterior = watershed_ABC_100k_obj$n2_posteriors$watershed_n2_posterior[is_noncoding_1,1],
              gam_n2_posterior = watershed_ABC_100k_obj$n2_posteriors$gam_n2_posterior[is_noncoding_1,1],
              rvot1 = real_valued_outliers_test1[is_noncoding_1,1],
              bot2 = binary_outliers_test2_alt[is_noncoding_2,1])
    bootstrap_delta_auprc <- boot(eval_data,statistic = stat_compute_delta_auprc, R = 5000)
    bootstrap_ci_delta_auprc <- boot.ci(bootstrap_delta_auprc)
    return (list(boot_out = bootstrap_delta_auprc, boot_ci = bootstrap_ci_delta_auprc))
  }else {
    noncoding_watershed <- matrix(watershed_ABC_100k_obj$n2_posteriors$watershed_n2_posterior[is_noncoding_1,])
    noncoding_gam <- matrix(watershed_ABC_100k_obj$n2_posteriors$gam_n2_posterior[is_noncoding_1,])
    noncoding_rvot1 <- matrix(real_valued_outliers_test1[is_noncoding_1,])
    noncoding_bot2 <- matrix(binary_outliers_test2_alt[is_noncoding_2,])
    auc_object_across_dimensions <- compute_auc_across_dimensions(1,
                                                              matrix(noncoding_watershed),
                                                              matrix(noncoding_gam), 
                                                              matrix(noncoding_rvot1), 
                                                              noncoding_bot2)
    return(auc_object_across_dimensions)
  }
}


make_panel_medz <- function(df) {
    # Set ordering of model types
	df$prediction_type=factor(df$prediction_type, levels=c("Watershed", "GAM"))

    # Make Precision-recall curve for ASE
    # This plot contains 3 curves (one for Watershed, RIVER, and GAM)
	outlier_type <- 'Median'
	outlier_name <- "Multiple Tissue Median Zscore"
  	plotter_ase <- ggplot(data=df[as.character(df$outlier_type)==outlier_type,], aes(x=recall, y=precision, group=prediction_type)) + geom_line(aes(linetype=prediction_type, colour=prediction_type)) + 
                labs(x="Recall", y="Precision", group="", linetype="", colour="") +
                scale_y_continuous(expand = c(0, 0), limits = c(0, 1)) + 
                scale_x_continuous(expand = c(0, 0), limits = c(0, 1)) + 
                theme(legend.position="bottom") +
                theme(panel.spacing = unit(2, "lines")) +
                scale_color_manual(values=c("#2C70B0","#FFBC7B")) +
                scale_linetype_manual(values=c("solid", "solid")) +
                gtex_v8_figure_theme() + 
                draw_label(outlier_name,x=.5,y=.95,size=8)
  	return (plotter_ase)
}

make_panel_4d <- function(df) {
    # Set ordering of model types
	df$prediction_type=factor(df$prediction_type, levels=c("Watershed", "GAM"))

    # Make Precision-recall curve for ASE
    # This plot contains 3 curves (one for Watershed, RIVER, and GAM)
	outlier_type <- 'Blood'
	outlier_name <- "Blood"
  	plotter_ase <- ggplot(data=df[as.character(df$outlier_type)==outlier_type,], aes(x=recall, y=precision, group=prediction_type)) + geom_line(aes(linetype=prediction_type, colour=prediction_type)) + 
                labs(x="Recall", y="Precision", group="", linetype="", colour="") +
                scale_y_continuous(expand = c(0, 0), limits = c(0, 1)) + 
                scale_x_continuous(expand = c(0, 0), limits = c(0, 1)) + 
                theme(legend.position="bottom") +
                theme(panel.spacing = unit(2, "lines")) +
                scale_color_manual(values=c("#2C70B0","#FFBC7B")) +
                scale_linetype_manual(values=c("solid", "solid")) +
                gtex_v8_figure_theme() + 
                draw_label(outlier_name,x=.5,y=.95,size=8)

    # Make Precision-recall curve for splicing
    # This plot contains 3 curves (one for Watershed, RIVER, and GAM)
	outlier_type <- "Fibroblast"
	outlier_name <- "Fibroblast"
  	plotter_splice <- ggplot(data=df[as.character(df$outlier_type)==outlier_type,], aes(x=recall, y=precision, group=prediction_type)) + geom_line(aes(linetype=prediction_type, colour=prediction_type)) + 
                labs(x="Recall", y="Precision",group="") +
                scale_y_continuous(expand = c(0, 0), limits = c(0, 1)) + 
                scale_x_continuous(expand = c(0, 0), limits = c(0, 1)) + 
                theme(legend.position="bottom") +
                theme(panel.spacing = unit(2, "lines")) +
                scale_color_manual(values=c("#2C70B0","#FFBC7B")) +
                scale_linetype_manual(values=c("solid", "solid")) +
                gtex_v8_figure_theme() +
                 draw_label(outlier_name,x=.5,y=.95,size=8)

    # Make Precision-recall curve for Expression
    # This plot contains 3 curves (one for Watershed, RIVER, and GAM)
	outlier_type <- "Muscle"
	outlier_name <- "Muscle"
  	plotter_expression <- ggplot(data=df[as.character(df$outlier_type)==outlier_type,], aes(x=recall, y=precision, group=prediction_type)) + geom_line(aes(linetype=prediction_type, colour=prediction_type)) + 
                labs(x="Recall", y="Precision", group="") +
                scale_y_continuous(expand = c(0, 0), limits = c(0, 1)) + 
                scale_x_continuous(expand = c(0, 0), limits = c(0, 1)) + 
                theme(legend.position="bottom") +
                theme(panel.spacing = unit(2, "lines")) +
                scale_color_manual(values=c("#2C70B0","#FFBC7B")) +
                scale_linetype_manual(values=c("solid", "solid")) +
                gtex_v8_figure_theme() + 
                draw_label(outlier_name,x=.5,y=.95,size=8)
  	# Make Precision-recall curve for Expression
    # This plot contains 3 curves (one for Watershed, RIVER, and GAM)
	outlier_type <- "MedZ"
	outlier_name <- "MedZ"
  	plotter_medZ <- ggplot(data=df[as.character(df$outlier_type)==outlier_type,], aes(x=recall, y=precision, group=prediction_type)) + geom_line(aes(linetype=prediction_type, colour=prediction_type)) + 
                labs(x="Recall", y="Precision", group="") +
                scale_y_continuous(expand = c(0, 0), limits = c(0, 1)) + 
                scale_x_continuous(expand = c(0, 0), limits = c(0, 1)) + 
                theme(legend.position="bottom") +
                theme(panel.spacing = unit(2, "lines")) +
                scale_color_manual(values=c("#2C70B0","#FFBC7B")) +
                scale_linetype_manual(values=c("solid", "solid")) +
                gtex_v8_figure_theme() + 
                draw_label(outlier_name,x=.5,y=.95,size=8)

    # Get legend from ASE Plot (Note: this is the same legend as found in the the Expression and Splicing plots)
    legend <- get_legend(plotter_ase + theme(legend.position="bottom"))

    # Combine PR-curves from Expression, ASE, and splicing into 1 plot via cowplot
    combined_plots <- plot_grid(plotter_expression + theme(legend.position="none"), 
                                plotter_ase + theme(legend.position="none"), 
                                plotter_splice+ theme(legend.position="none"), 
                                plotter_medZ+ theme(legend.position="none"), 
                                rel_widths=c(1,1,1,1), nrow=2)

    # Merge combined plot with single legend via cowplot
	combined <- ggdraw() + draw_plot(combined_plots,0,.07,1,.9) + draw_plot(legend,.38,-.40,1,1)
	
	return(combined)
}

```
## Figure 3.1 Watershed-SV outperforms WGS baseline method. 
```{r}
watershed_ABC_100k_obj_Cells_Cultured_Fibroblasts_file  <- '/oak/stanford/groups/smontgom/tannerj/UDN_LRS_watershed/Watershed-SV/new-protein-lincRNA-100k-rare-position-aware-2.0/eval_watarshed_ABC_regional_Cells_Cultured_fibroblast_SV_RIVER_100000_number_of_dimensions_1_pthresh_0.000063_pval_frac_0.176_separate_trainP_2.56__evaluation_object.rds'
data_Cells_Cultured_Fibroblasts_file <- '/oak/stanford/groups/smontgom/tannerj/UDN_LRS_watershed/Watershed-SV/new-protein-lincRNA-100k-rare-position-aware-2.0/eval_watarshed_ABC_regional.Cells_Cultured_fibroblasts.z4.tsv'
watershed_ABC_100k_obj_Muscle_Skeletal_file <- '/oak/stanford/groups/smontgom/tannerj/UDN_LRS_watershed/Watershed-SV/new-protein-lincRNA-100k-rare-position-aware-2.0/eval_watarshed_ABC_regional_Muscle_Skeletal_SV_RIVER_100000_number_of_dimensions_1_pthresh_0.000063_pval_frac_0.151_separate_trainP_2.56__evaluation_object.rds'
data_Muscle_Skeletal_file <- '/oak/stanford/groups/smontgom/tannerj/UDN_LRS_watershed/Watershed-SV/new-protein-lincRNA-100k-rare-position-aware-2.0/eval_watarshed_ABC_regional.Muscle_Skeletal.z4.tsv'
watershed_ABC_100k_obj_Whole_Blood_file <- '/oak/stanford/groups/smontgom/tannerj/UDN_LRS_watershed/Watershed-SV/new-protein-lincRNA-100k-rare-position-aware-2.0/eval_watarshed_ABC_regional_Whole_Blood_SV_RIVER_100000_number_of_dimensions_1_pthresh_0.000063_pval_frac_0.138_separate_trainP_2.56__evaluation_object.rds'
data_Whole_Blood_file <- '/oak/stanford/groups/smontgom/tannerj/UDN_LRS_watershed/Watershed-SV/new-protein-lincRNA-100k-rare-position-aware-2.0/eval_watarshed_ABC_regional.Whole_Blood.z4.tsv'
watershed_ABC_100k_obj_Median_file  <- '/oak/stanford/groups/smontgom/tannerj/UDN_LRS_watershed/Watershed-SV/new-protein-lincRNA-100k-rare-position-aware-2.0/eval_watarshed_ABC_regional_MedZ_SV_RIVER_100000_number_of_dimensions_1_pthresh_0.0027_pval_frac_0.16_separate_trainP_2.56__evaluation_object.rds'
data_Median_file <- '/oak/stanford/groups/smontgom/tannerj/UDN_LRS_watershed/Watershed-SV/new-protein-lincRNA-100k-rare-position-aware-2.0/eval_watarshed_ABC_regional.medZ.z3.tsv'
watershed_ABC_100k_obj_Median_baseline_file  <- '/oak/stanford/groups/smontgom/tannerj/UDN_LRS_watershed/Watershed-SV/new-protein-lincRNA-100k-rare-baseline/eval_watershed_baseline_SV_RIVER_100000_number_of_dimensions_1_pthresh_0.0027_pval_frac_0.1618_shared_trainP_2.56__evaluation_object.rds'
data_Median_baseline_file <- '/oak/stanford/groups/smontgom/tannerj/UDN_LRS_watershed/Watershed-SV/new-protein-lincRNA-100k-rare-baseline/eval_watarshed_ABC_regional.medz.z3.tsv'
watershed_ABC_10k_obj_Median_file  <- '/oak/stanford/groups/smontgom/tannerj/UDN_LRS_watershed/Watershed-SV/new-protein-lincRNA-10k-final-version/eval_10k_annotations_SV_RIVER_10000_number_of_dimensions_1_pthresh_0.0027_pval_frac_0.31_separate_trainP_2.56__evaluation_object.rds'
data_Median_10k_file <- '/oak/stanford/groups/smontgom/tannerj/UDN_LRS_watershed/Watershed-SV/new-protein-lincRNA-10k-final-version/eval_10k_annotations.medZ.tsv'
```
# coding
```{r}
fibroblast_100k_bootstrap_test <- re_calculate_precision_recall(watershed_ABC_100k_obj_Cells_Cultured_Fibroblasts_file,data_Cells_Cultured_Fibroblasts_file,1,0.176,0.01,compute_ci = TRUE)
blood_100k_bootstrap_test <- re_calculate_precision_recall(watershed_ABC_100k_obj_Whole_Blood_file,data_Whole_Blood_file,1,0.176,0.01,compute_ci = TRUE)
muscle_100k_bootstrap_test <- re_calculate_precision_recall(watershed_ABC_100k_obj_Muscle_Skeletal_file,data_Muscle_Skeletal_file,1,0.176,0.01,compute_ci = TRUE)

```
```{r}
median_100k_bootstrap_test <- re_calculate_precision_recall(watershed_ABC_100k_obj_Median_file,data_Median_file,1,0.176,0.05,compute_ci = F)
```
```{r}
median_100k_baseline_bootstrap_test <- re_calculate_precision_recall(watershed_ABC_100k_obj_Median_baseline_file,data_Median_baseline_file,1,0.1618,0.045,compute_ci = F)
```

```{r}
median_10k_bootstrap_test <- re_calculate_precision_recall(watershed_ABC_10k_obj_Median_file,data_Median_10k_file,1,0.31,0.05,compute_ci = FALSE)
```
```{r}
median_100k_bootstrap_test$boot_out
```

```{r}
median_100k_bootstrap_test$boot_ci
```

```{r}
median_10k_bootstrap_test$boot_out
median_10k_bootstrap_test$boot_ci

```
```{r}
make_panel_medz_alt <- function(df) {
    # Set ordering of model types
	df$prediction_type=factor(df$prediction_type, levels=c("Watershed-10kb", 
	                                                       "Watershed-100kb", 
	                                                       "WGS-only-10kb",
	                                                       "WGS-only-100kb"))

    # Make Precision-recall curve for ASE
    # This plot contains 3 curves (one for Watershed, RIVER, and GAM)
	outlier_type <- 'Median'
	outlier_name <- "Multiple Tissue Median Zscore"
  	plotter_ase <- ggplot(data=df[as.character(df$outlier_type)==outlier_type,], aes(x=recall, y=precision, group=prediction_type)) + geom_line(aes(linetype=prediction_type, colour=prediction_type)) + 
                labs(x="Recall", y="Precision", group="", linetype="", colour="") +
                scale_y_continuous(expand = c(0, 0), limits = c(0, 1)) + 
                scale_x_continuous(expand = c(0, 0), limits = c(0, 1)) + 
                theme(legend.position="bottom") +
                theme(panel.spacing = unit(2, "lines")) +
                scale_color_manual(values=c("#2C70B0","#2C70B0","#FFBC7B","#FFBC7B")) +
                scale_linetype_manual(values=c("solid", "dashed","solid","dashed")) +
                gtex_v8_figure_theme() + 
                ggtitle(outlier_name,)
  	return (plotter_ase)
}

```

```{r}

tissue_dim <- 1
medz_watershed_100kb <- data_frame(recall = median_100k_bootstrap_test[[tissue_dim]]$evaROC$watershed_recall,
           precision = median_100k_bootstrap_test[[tissue_dim]]$evaROC$watershed_precision,
           prediction_type = 'Watershed-100kb',outlier_type = 'Median')
medz_gam_100kb <- data_frame(recall = median_100k_bootstrap_test[[tissue_dim]]$evaROC$GAM_recall,
           precision = median_100k_bootstrap_test[[tissue_dim]]$evaROC$GAM_precision,
           prediction_type = 'WGS-only-100kb',outlier_type = 'Median')
medz_watershed_10kb <- data_frame(recall = median_10k_bootstrap_test[[tissue_dim]]$evaROC$watershed_recall,
           precision = median_10k_bootstrap_test[[tissue_dim]]$evaROC$watershed_precision,
           prediction_type = 'Watershed-10kb',outlier_type = 'Median')
medz_gam_10kb <- data_frame(recall = median_10k_bootstrap_test[[tissue_dim]]$evaROC$GAM_recall,
           precision = median_10k_bootstrap_test[[tissue_dim]]$evaROC$GAM_precision,
           prediction_type = 'WGS-only-10kb',outlier_type = 'Median')
total_medz_pr_data <- rbind(medz_watershed_100kb,medz_gam_100kb,medz_watershed_10kb,medz_gam_10kb)
p <- make_panel_medz_alt(total_medz_pr_data) +theme(
        axis.text.x  = element_text(size = 11),
        axis.text.y  = element_text(size = 11),
        plot.title = element_text(size = 12),
        plot.margin = margin(2, 2, 2, 2, "cm")
  )

ggsave('Figure3/MedianZ_10kb_100kb.svg',p,height = 6,width = 6) 
```
```{r}
p
```
```{r}
tissue_dim <- 1
medz_watershed <- data_frame(recall = median_100k_baseline_bootstrap_test[[tissue_dim]]$evaROC$watershed_recall,
           precision = median_100k_baseline_bootstrap_test[[tissue_dim]]$evaROC$watershed_precision,
           prediction_type = 'Watershed',outlier_type = 'Median')
medz_gam <- data_frame(recall = median_100k_baseline_bootstrap_test[[tissue_dim]]$evaROC$GAM_recall,
           precision = median_100k_baseline_bootstrap_test[[tissue_dim]]$evaROC$GAM_precision,
           prediction_type = 'GAM',outlier_type = 'Median')
total_medz_pr_data <- rbind(medz_watershed,medz_gam)
p <- make_panel_medz(total_medz_pr_data) 
```
```{r}
p
```

```{r}
watershed_ABC_100k_obj_Cells_Cultured_Fibroblasts <- re_calculate_precision_recall(watershed_ABC_100k_obj_Cells_Cultured_Fibroblasts_file,data_Cells_Cultured_Fibroblasts_file,1,0.176,0.05)
watershed_ABC_100k_obj_Muscle_Skeletal <- re_calculate_precision_recall(watershed_ABC_100k_obj_Muscle_Skeletal_file,data_Muscle_Skeletal_file,1,0.152,0.05)
watershed_ABC_100k_obj_Whole_Blood <- re_calculate_precision_recall(watershed_ABC_100k_obj_Whole_Blood_file,data_Whole_Blood_file,1,0.138,0.05)
watershed_ABC_100k_obj_MedZ <- re_calculate_precision_recall(watershed_ABC_100k_obj_Median_file,data_Median_file,1,0.176,0.05)


blood_watershed <- data_frame(recall = watershed_ABC_100k_obj_Whole_Blood[[1]]$evaROC$watershed_recall,
           precision = watershed_ABC_100k_obj_Whole_Blood[[1]]$evaROC$watershed_precision,
           prediction_type = 'Watershed',outlier_type = 'Blood')
fibroblast_watershed <- data_frame(recall = watershed_ABC_100k_obj_Cells_Cultured_Fibroblasts[[1]]$evaROC$watershed_recall,
           precision = watershed_ABC_100k_obj_Cells_Cultured_Fibroblasts[[1]]$evaROC$watershed_precision,
           prediction_type = 'Watershed',outlier_type = 'Fibroblast')
muscle_watershed <- data_frame(recall = watershed_ABC_100k_obj_Muscle_Skeletal[[1]]$evaROC$watershed_recall,
           precision = watershed_ABC_100k_obj_Muscle_Skeletal[[1]]$evaROC$watershed_precision,
           prediction_type = 'Watershed',outlier_type = 'Muscle')
medz_watershed <- data_frame(recall = watershed_ABC_100k_obj_MedZ[[1]]$evaROC$watershed_recall,
           precision = watershed_ABC_100k_obj_MedZ[[1]]$evaROC$watershed_precision,
           prediction_type = 'Watershed',outlier_type = 'MedZ')
blood_gam <- data_frame(recall = watershed_ABC_100k_obj_Whole_Blood[[1]]$evaROC$GAM_recall,
           precision = watershed_ABC_100k_obj_Whole_Blood[[1]]$evaROC$GAM_precision,
           prediction_type = 'GAM',outlier_type = 'Blood')
fibroblast_gam <- data_frame(recall = watershed_ABC_100k_obj_Cells_Cultured_Fibroblasts[[1]]$evaROC$GAM_recall,
           precision = watershed_ABC_100k_obj_Cells_Cultured_Fibroblasts[[1]]$evaROC$GAM_precision,
           prediction_type = 'GAM',outlier_type = 'Fibroblast')
muscle_gam <- data_frame(recall = watershed_ABC_100k_obj_Muscle_Skeletal[[1]]$evaROC$GAM_recall,
           precision = watershed_ABC_100k_obj_Muscle_Skeletal[[1]]$evaROC$GAM_precision,
           prediction_type = 'GAM',outlier_type = 'Muscle')
medz_gam <- data_frame(recall = watershed_ABC_100k_obj_MedZ[[1]]$evaROC$GAM_recall,
           precision = watershed_ABC_100k_obj_MedZ[[1]]$evaROC$GAM_precision,
           prediction_type = 'GAM',outlier_type = 'MedZ')
total_pr_data <- rbind(blood_watershed,fibroblast_watershed,muscle_watershed,medz_watershed,blood_gam,fibroblast_gam,muscle_gam,medz_gam)

fig_3_supplementary <- make_panel_4d(total_pr_data) + theme(plot.margin = margin(2,2,2,2,"cm"))

ggsave('Figure3/Supplementary_6a.tissue_watershed_medz_watershed.pdf',fig_3_supplementary,height = 6,width = 6)
```

```{r}
watershed_ABC_100k_obj_Muscle_Skeletal[[1]]$evaROC$watershed_pr_auc
```

# noncoding
```{r}
watershed_ABC_100k_obj_Cells_Cultured_Fibroblasts_noncoding <- re_calculate_precision_recall_noncoding(watershed_ABC_100k_obj_Cells_Cultured_Fibroblasts_file,data_Cells_Cultured_Fibroblasts_file,1,0.2,0.1,compute_ci = FALSE)
watershed_ABC_100k_obj_Muscle_Skeletal_noncoding <- re_calculate_precision_recall_noncoding(watershed_ABC_100k_obj_Muscle_Skeletal_file,data_Muscle_Skeletal_file,1,0.2,0.05,compute_ci = FALSE)
watershed_ABC_100k_obj_Whole_Blood_noncoding <- re_calculate_precision_recall_noncoding(watershed_ABC_100k_obj_Whole_Blood_file,data_Whole_Blood_file,1,0.2,0.1,compute_ci = FALSE)
watershed_ABC_100k_obj_MedZ_noncoding <- re_calculate_precision_recall_noncoding(watershed_ABC_100k_obj_Median_file,data_Median_file,1,0.2,0.1,compute_ci = FALSE)


blood_watershed_noncoding <- data_frame(recall = watershed_ABC_100k_obj_Whole_Blood_noncoding[[1]]$evaROC$watershed_recall,
           precision = watershed_ABC_100k_obj_Whole_Blood_noncoding[[1]]$evaROC$watershed_precision,
           prediction_type = 'Watershed',outlier_type = 'Blood')
fibroblast_watershed_noncoding <- data_frame(recall = watershed_ABC_100k_obj_Cells_Cultured_Fibroblasts_noncoding[[1]]$evaROC$watershed_recall,
           precision = watershed_ABC_100k_obj_Cells_Cultured_Fibroblasts_noncoding[[1]]$evaROC$watershed_precision,
           prediction_type = 'Watershed',outlier_type = 'Fibroblast')
muscle_watershed_noncoding <- data_frame(recall = watershed_ABC_100k_obj_Muscle_Skeletal_noncoding[[1]]$evaROC$watershed_recall,
           precision = watershed_ABC_100k_obj_Muscle_Skeletal_noncoding[[1]]$evaROC$watershed_precision,
           prediction_type = 'Watershed',outlier_type = 'Muscle')
medz_watershed_noncoding <- data_frame(recall = watershed_ABC_100k_obj_MedZ_noncoding[[1]]$evaROC$watershed_recall,
           precision = watershed_ABC_100k_obj_MedZ_noncoding[[1]]$evaROC$watershed_precision,
           prediction_type = 'Watershed',outlier_type = 'MedZ')
blood_gam_noncoding <- data_frame(recall = watershed_ABC_100k_obj_Whole_Blood_noncoding[[1]]$evaROC$GAM_recall,
           precision = watershed_ABC_100k_obj_Whole_Blood_noncoding[[1]]$evaROC$GAM_precision,
           prediction_type = 'GAM',outlier_type = 'Blood')
fibroblast_gam_noncoding <- data_frame(recall = watershed_ABC_100k_obj_Cells_Cultured_Fibroblasts_noncoding[[1]]$evaROC$GAM_recall,
           precision = watershed_ABC_100k_obj_Cells_Cultured_Fibroblasts_noncoding[[1]]$evaROC$GAM_precision,
           prediction_type = 'GAM',outlier_type = 'Fibroblast')
muscle_gam_noncoding <- data_frame(recall = watershed_ABC_100k_obj_Muscle_Skeletal_noncoding[[1]]$evaROC$GAM_recall,
           precision = watershed_ABC_100k_obj_Muscle_Skeletal_noncoding[[1]]$evaROC$GAM_precision,
           prediction_type = 'GAM',outlier_type = 'Muscle')
medz_gam_noncoding <- data_frame(recall = watershed_ABC_100k_obj_MedZ_noncoding[[1]]$evaROC$GAM_recall,
           precision = watershed_ABC_100k_obj_MedZ_noncoding[[1]]$evaROC$GAM_precision,
           prediction_type = 'GAM',outlier_type = 'MedZ')
total_pr_data_noncoding <- rbind(fibroblast_watershed_noncoding,blood_watershed_noncoding,muscle_watershed_noncoding,medz_watershed_noncoding,
                                 fibroblast_gam_noncoding,blood_gam_noncoding,muscle_gam_noncoding,medz_gam_noncoding)

p_noncoding <- make_panel_4d(total_pr_data_noncoding) + labs('Noncoding performance in different outlier signals') + theme(plot.margin = margin(2,2,2,2,"cm"))
ggsave('Figure3/Supplementary_6b.noncoding.tissue_watershed_medz_watershed.pdf',p_noncoding,height = 6,width = 6)

```
```{r}
p_noncoding
```

```{r}
median_10k_obj <- readRDS('/oak/stanford/groups/smontgom/tannerj/UDN_LRS_watershed/Watershed-SV/new-protein-lincRNA-10k-final-version/eval_10k_annotations_SV_RIVER_10000_number_of_dimensions_1_pthresh_0.0027_pval_frac_0.31_separate_trainP_2.56__evaluation_object.rds')
annotations_train <- read_delim('/oak/stanford/groups/smontgom/tannerj/UDN_LRS_watershed/Watershed-SV/new-protein-lincRNA-10k-final-version/eval_10k_annotations.medZ.tsv',delim = '\t')
feature_categories <-read.csv('./collapse_annotation_instructions.csv')
feature_names_pre_filter <- annotations_train %>% colnames()
feature_names <- feature_names_pre_filter[!feature_names_pre_filter %in% c("SubjectID","GeneName","TE_pvalues_","N2pair")]
feature_categories_nona <- feature_categories %>% filter(Annotation.Type!='') 
```

```{r}
theta_medz_df_watershed <- data.frame(feature = feature_names, 
                       beta = median_10k_obj$model_params$theta,
                       tissue = 'Median',model='Watershed')
theta_medz_df_gam <- data.frame(feature = feature_names, 
                       beta = median_10k_obj$gam_model_params$gam_parameters$theta,
                       tissue = 'Median',model='supervised baseline')
theta_medz_df <- rbind(theta_medz_df_watershed,theta_medz_df_gam)

total_theta_df_w_categories <- theta_medz_df %>% left_join(feature_categories_nona,by = c("feature" = "Feature")) %>% dplyr::select(-region.revision)
```
```{r}
all_features <- total_theta_df_w_categories %>% 
  mutate(tissue = as.factor(tissue), Annotation.Type = as.factor(Annotation.Type), region = as.factor(region)) %>%
  filter((tissue=='Median') & (model=='Watershed')) %>% 
  group_by(Annotation.Type) 
effect_plot_median_10k <- all_features %>% 
  ggplot(aes(x = beta, y = reorder(feature,-beta), fill = Annotation.Type)) + 
  geom_bar(stat = 'identity', position = "dodge", ) + ylab( label = 'annotations' ) + xlab( label = 'effect size' ) + 
  facet_grid(rows = vars(Annotation.Type),cols = vars(model), scales="free", space="free_y") 

total_theta_df_w_categories$model <- factor(total_theta_df_w_categories$model,levels = c('Watershed','supervised baseline'))
effect_plot_median_10k

effect_plot_median_10k %>% ggsave(file="Median_10kb_baseline_gam_watershed_effectsize_final.pdf", width=10, height=20, dpi=300)
```

```{r}
median_100k_obj <- readRDS('/oak/stanford/groups/smontgom/tannerj/UDN_LRS_watershed/Watershed-SV/new-protein-lincRNA-100k-rare-position-aware-2.0/eval_watarshed_ABC_regional_MedZ_SV_RIVER_100000_number_of_dimensions_1_pthresh_0.0027_pval_frac_0.163_separate_trainP_2.56__evaluation_object.rds')
feature_categories <-read.csv('./collapse_annotation_instructions.csv')
feature_names <- as.data.frame(median_100k_obj$feat_test) %>% colnames()
feature_categories_nona <- feature_categories %>% filter(Annotation.Type!='') 

theta_medz_df_watershed <- data.frame(feature = feature_names, 
                       beta = median_100k_obj$model_params$theta,
                       tissue = 'Median',model='Watershed')
theta_medz_df_gam <- data.frame(feature = feature_names, 
                       beta = median_100k_obj$gam_model_params$gam_parameters$theta,
                       tissue = 'Median',model='supervised baseline')
theta_medz_df <- rbind(theta_medz_df_watershed,theta_medz_df_gam)

total_theta_df_w_categories <- theta_medz_df %>% left_join(feature_categories_nona,by = c("feature" = "Feature")) %>% dplyr::select(-region.revision)
```
```{r}
total_theta_df_w_categories$region[total_theta_df_w_categories$region == 'None'] <- 'Nonregion-specific'
top_positive_features <- total_theta_df_w_categories %>% 
  mutate(tissue = as.factor(tissue), Annotation.Type = as.factor(Annotation.Type), region = as.factor(region)) %>%
  filter((tissue=='Median') & (model=='Watershed')) %>% 
  group_by(Annotation.Type,region) 

#effect_plot_median_100k <- rbind(top_positive_features,top_negative_features) %>% 
effect_plot_median_100k <- top_positive_features %>% 
  ggplot(aes(x = beta, y = reorder(feature,-beta), fill = Annotation.Type)) + 
  geom_bar(stat = 'identity', position = "dodge", ) + ylab( label = 'annotations' ) + xlab( label = 'effect size' ) + 
  facet_grid(rows = vars(factor(region,levels = c('Nonregion-specific','tss flanking','Gene Body','tes flanking'))), scales="free", space="free_y") + theme(plot.title=element_text(face="bold"))

total_theta_df_w_categories$model <- factor(total_theta_df_w_categories$model,levels = c('Watershed','supervised baseline'))
effect_plot_median_100k

effect_plot_median_100k %>% ggsave(file="Median_100kb_region_gam_watershed_effectsize_reduced_final.pdf", width=15, height=20, dpi=300)
```

```{r}
top_positive_features
```
```{r}
median_100k_obj$model_params$phi
```

```{r}
median_100k_obj <- readRDS('/oak/stanford/groups/smontgom/tannerj/UDN_LRS_watershed/Watershed-SV/new-protein-lincRNA-10k-final-version/eval_10k_annotations_SV_RIVER_10000_number_of_dimensions_1_pthresh_0.0027_pval_frac_0.31_separate_trainP_2.56__evaluation_object.rds')
feature_categories <-read.csv('./collapse_annotation_instructions.csv')
feature_names <- as.data.frame(median_100k_obj$feat_test) %>% colnames()
feature_categories_nona <- feature_categories %>% filter(Annotation.Type!='') 

theta_medz_df_watershed <- data.frame(feature = feature_names, 
                       beta = median_100k_obj$model_params$theta,
                       tissue = 'Median',model='Watershed')
theta_medz_df_gam <- data.frame(feature = feature_names, 
                       beta = median_100k_obj$gam_model_params$gam_parameters$theta,
                       tissue = 'Median',model='supervised baseline')
theta_medz_df <- rbind(theta_medz_df_watershed,theta_medz_df_gam)

total_theta_df_w_categories <- theta_medz_df %>% left_join(feature_categories_nona,by = c("feature" = "Feature")) %>% dplyr::select(-region.revision)
```
```{r}
top_positive_features <- total_theta_df_w_categories %>% 
  mutate(tissue = as.factor(tissue), Annotation.Type = as.factor(Annotation.Type), region = as.factor(region)) %>%
  filter((tissue=='Median') & (model=='Watershed')) %>% 
  group_by(Annotation.Type) %>% 
  top_n(n = 1, beta) 
top_negative_features <- total_theta_df_w_categories %>% 
  mutate(tissue = as.factor(tissue), Annotation.Type = as.factor(Annotation.Type), region = as.factor(region)) %>%
  filter((tissue=='Median') & (model=='Watershed')) %>% 
  group_by(Annotation.Type) %>% 
  top_n(n = -1, beta) 
#effect_plot_median_10k <- rbind(top_positive_features,top_negative_features) %>% 
effect_plot_median_10k <- top_positive_features %>% 
  ggplot(aes(x = beta, y = reorder(feature,-beta), fill = Annotation.Type)) + 
  geom_bar(stat = 'identity', position = "dodge", ) + ylab( label = 'annotations' ) + xlab( label = 'effect size' )# + 
  #facet_grid(rows = vars(Annotation.Type),cols = vars(model), scales="free", space="free_y") 

total_theta_df_w_categories$model <- factor(total_theta_df_w_categories$model,levels = c('Watershed','supervised baseline'))
effect_plot_median_10k

effect_plot_median_10k %>% ggsave(file="Figure3/Median_10kb_baseline_gam_watershed_effectsize_reduced.svg", width=15, height=10, dpi=300)
```
# effect size from median z-score 100kb model. 
```{r}
median_100k_obj <- readRDS('/oak/stanford/groups/smontgom/tannerj/UDN_LRS_watershed/Watershed-SV/new-protein-lincRNA-100k-rare-position-aware-2.0/eval_watarshed_ABC_regional_MedZ_SV_RIVER_100000_number_of_dimensions_1_pthresh_0.0027_pval_frac_0.163_separate_trainP_2.56__evaluation_object.rds')
feature_categories <-read.csv('./collapse_annotation_instructions.csv')
feature_names <- as.data.frame(median_100k_obj$feat_test) %>% colnames()
feature_categories_nona <- feature_categories %>% filter(Annotation.Type!='') 

theta_medz_df_watershed <- data.frame(feature = feature_names, 
                       beta = median_100k_obj$model_params$theta,
                       tissue = 'Median',model='Watershed')
theta_medz_df_gam <- data.frame(feature = feature_names, 
                       beta = median_100k_obj$gam_model_params$gam_parameters$theta,
                       tissue = 'Median',model='supervised baseline')
theta_medz_df <- rbind(theta_medz_df_watershed,theta_medz_df_gam)

total_theta_df_w_categories <- theta_medz_df %>% left_join(feature_categories_nona,by = c("feature" = "Feature")) %>% dplyr::select(-region.revision)
```
```{r}
total_theta_df_w_categories$region[total_theta_df_w_categories$region == 'None'] <- 'Nonregion-specific'
top_positive_features <- total_theta_df_w_categories %>% 
  mutate(tissue = as.factor(tissue), Annotation.Type = as.factor(Annotation.Type), region = as.factor(region)) %>%
  filter((tissue=='Median') & (model=='Watershed')) %>% 
  group_by(Annotation.Type,region) %>% 
  top_n(n = 1, beta) 
top_negative_features <- total_theta_df_w_categories %>% 
  mutate(tissue = as.factor(tissue), Annotation.Type = as.factor(Annotation.Type), region = as.factor(region)) %>%
  filter((tissue=='Median') & (model=='Watershed')) %>% 
  group_by(Annotation.Type) %>% 
  top_n(n = -1, beta) 
#effect_plot_median_100k <- rbind(top_positive_features,top_negative_features) %>% 
effect_plot_median_100k <- top_positive_features %>% 
  ggplot(aes(x = beta, y = reorder(feature,-beta), fill = Annotation.Type)) + 
  geom_bar(stat = 'identity', position = "dodge", ) + ylab( label = 'annotations' ) + xlab( label = 'effect size' ) + 
  facet_grid(rows = vars(factor(region,levels = c('Nonregion-specific','tss flanking','Gene Body','tes flanking'))), scales="free", space="free_y") + theme(plot.title=element_text(face="bold"))

total_theta_df_w_categories$model <- factor(total_theta_df_w_categories$model,levels = c('Watershed','supervised baseline'))
effect_plot_median_100k

effect_plot_median_100k %>% ggsave(file="Figure3/Median_100kb_region_gam_watershed_effectsize_reduced.svg", width=15, height=20, dpi=300)
```
```{r}
top_positive_features
```

# plot feature contributions
```{r}
total_theta_df_w_categories$model <- factor(total_theta_df_w_categories$model,levels = c('Watershed','supervised baseline'))
```

```{r}
effect_plot <- total_theta_df_w_categories %>% mutate(tissue = as.factor(tissue),Annotation.Type = as.factor(Annotation.Type), region = as.factor(region)) %>% filter(tissue=='Median') %>% ggplot(aes(x = beta, y = reorder(feature,-beta), fill = Annotation.Type)) + 
  geom_bar(stat = 'identity', position = "dodge", ) + 
  facet_grid(cols = vars(tissue)) 
effect_plot %>% ggsave(file="Baseline_watershed_SV_effectsize.w_vep3.10kb.medZ.png", width=6, height=15, dpi=300)
```

```{r}
watershed_ABC_100k_obj_Cells_Cultured_Fibroblasts_noncoding[[1]]$evaROC$GAM_pr_auc
```

```{r}
watershed_ABC_100k_obj <-readRDS('/oak/stanford/groups/smontgom/tannerj/UDN_LRS_watershed/Watershed-SV/new-protein-lincRNA-10k-final-version/eval_10k_annotations_SV_RIVER_10000_number_of_dimensions_1_pthresh_0.0027_pval_frac_0.31_separate_trainP_2.56__evaluation_object.rds')
input_file <- '/oak/stanford/groups/smontgom/tannerj/UDN_LRS_watershed/Watershed-SV/new-protein-lincRNA-10k-final-version/eval_10k_annotations.medZ.tsv'
data_input <- load_watershed_data(input_file, 1, 0.32, 0.05)
# Parse data_input for evaluation-related relevant fields
feat_all <- data_input$feat
discrete_outliers_all <- data_input$outliers_discrete
binary_outliers_all <- data_input$outliers_binary
fraction_binary_outliers_all <- data_input$fraction_outliers_binary
N2_pairs <- data_input$N2_pairs

# Extraction of test data for Watershed
feat_test <- rbind(feat_all[!is.na(N2_pairs),][seq(from=1,to=sum(!is.na(N2_pairs)),by=2),], feat_all[!is.na(N2_pairs),][seq(from=2,to=sum(!is.na(N2_pairs)),by=2),])
discrete_outliers_test1 <- as.matrix(c(discrete_outliers_all[!is.na(N2_pairs)][seq(from=1,to=sum(!is.na(N2_pairs)),by=2)], discrete_outliers_all[!is.na(N2_pairs)][seq(from=2,to=sum(!is.na(N2_pairs)),by=2)]))

discrete_outliers_test2 <- as.matrix(c(discrete_outliers_all[!is.na(N2_pairs)][seq(from=2,to=sum(!is.na(N2_pairs)),by=2)], discrete_outliers_all[!is.na(N2_pairs)][seq(from=1,to=sum(!is.na(N2_pairs)),by=2)]))

binary_outliers_test1 <- as.matrix(c(binary_outliers_all[!is.na(N2_pairs)][seq(from=1,to=sum(!is.na(N2_pairs)),by=2)], binary_outliers_all[!is.na(N2_pairs)][seq(from=2,to=sum(!is.na(N2_pairs)),by=2)]))

binary_outliers_test2 <- as.matrix(c(fraction_binary_outliers_all[!is.na(N2_pairs)][seq(from=2,to=sum(!is.na(N2_pairs)),by=2)], fraction_binary_outliers_all[!is.na(N2_pairs)][seq(from=1,to=sum(!is.na(N2_pairs)),by=2)]))
		# Absolute pvalues from test prediction data set (to be used for RNA-only analysis)
real_valued_outliers_test1 <- -log10(abs(as.matrix(c(data_input$outlier_pvalues[!is.na(N2_pairs)][seq(from=1,to=sum(!is.na(N2_pairs)),by=2)], data_input$outlier_pvalues[!is.na(N2_pairs)][seq(from=2,to=sum(!is.na(N2_pairs)),by=2)]))) + 1e-7)

binary_outliers_test2_alt <- as.matrix(c(binary_outliers_all[!is.na(N2_pairs)][seq(from=2,to=sum(!is.na(N2_pairs)),by=2)], binary_outliers_all[!is.na(N2_pairs)][seq(from=1,to=sum(!is.na(N2_pairs)),by=2)]))

watershed <- watershed_ABC_100k_obj$n2_posteriors$watershed_n2_posterior
gam <- watershed_ABC_100k_obj$n2_posteriors$gam_n2_posterior
rvot1 <- real_valued_outliers_test1
bot2 <- binary_outliers_test2_alt
auc_object_across_dimensions <- compute_auc_across_dimensions(1,
                                                              watershed,
                                                              gam, 
                                                              rvot1, 
                                                              bot2)
```
```{r}
tissue_dim <- 1
medz_watershed <- data_frame(recall = auc_object_across_dimensions[[tissue_dim]]$evaROC$watershed_recall,
           precision = auc_object_across_dimensions[[tissue_dim]]$evaROC$watershed_precision,
           prediction_type = 'Watershed',outlier_type = 'Median')
medz_gam <- data_frame(recall = auc_object_across_dimensions[[tissue_dim]]$evaROC$GAM_recall,
           precision = auc_object_across_dimensions[[tissue_dim]]$evaROC$GAM_precision,
           prediction_type = 'GAM',outlier_type = 'Median')
total_medz_pr_data <- rbind(medz_watershed,medz_gam)
make_panel_medz(total_medz_pr_data)
```
# extract 10k model predictions and overlap with constraints and disease variant, gene set. 
```{r}
Watershed_10k_medz_pred_file <- '/oak/stanford/groups/smontgom/tannerj/UDN_LRS_watershed/Watershed-SV/new-protein-lincRNA-10k-final-version/MedZ_prediction_SV_RIVER_10000_number_of_dimensions_1_pthresh_0.0027_trainP_2.56__posterior_probability.txt'
Watershed_10k_CMG_prediction_file <- '/oak/stanford/groups/smontgom/tannerj/UDN_LRS_watershed/Watershed-SV/MuscularDystrophy-protein-lincRNA-100k-all-position-aware-3.0/MedZ_model_prediction_z3_SV_RIVER_100000_number_of_dimensions_1_pthresh_0.0027_trainP_2.56__posterior_probability.txt'
Watershed_10k_medz_pred <- read_delim(Watershed_10k_medz_pred_file,delim = '\t')
Watershed_10k_cmg_medz_pred <- read_delim(Watershed_10k_CMG_prediction_file,delim = '\t')

#gnomAD_v4_constraints_file <- '/oak/stanford/groups/smontgom/tannerj/UDN_LRS_watershed/Watershed-SV/input/gnomad.v4.0.constraint_metrics.tsv'
#gnomAD_v4_constraints <- read_delim(gnomAD_v4_constraints_file,delim = '\t')
Watershed_10k_medz_pred <- Watershed_10k_medz_pred %>% separate(sample_names, into = c("SubjectID","GeneName","SV"), sep = ":", extra = "merge")
Watershed_10k_cmg_medz_pred <- Watershed_10k_cmg_medz_pred %>% separate(sample_names, into = c("SubjectID","GeneName","SV"), sep = ":", extra = "merge")

```

```{r}
ensembl_ids <- unique(as.vector(Watershed_10k_medz_pred$GeneName))
symbols <- select(org.Hs.eg.db, keys=ensembl_ids, columns="SYMBOL", keytype="ENSEMBL")
ensembl_ids_cmg <- unique(as.vector(Watershed_10k_cmg_medz_pred$GeneName))
symbols_cmg <- select(org.Hs.eg.db, keys=ensembl_ids_cmg, columns="SYMBOL", keytype="ENSEMBL")


```
```{r}
symbols_cmg %>% filter(is.na(SYMBOL))
```

```{r}
corrected_symbols <- symbols %>% filter(!(SYMBOL %in% c("IPO11-LRRC70","ZNF765-ZNF761","GRIK1-AS2")))
corrected_symbols_cmg <- symbols_cmg %>% filter(!(SYMBOL %in% c("LINC00673")))

Watershed_10k_medz_pred_symbol <- Watershed_10k_medz_pred %>% left_join(corrected_symbols,by = c("GeneName" = "ENSEMBL")) %>% dplyr::rename( gene = SYMBOL)
Watershed_10k_cmg_medz_pred_symbol <- Watershed_10k_cmg_medz_pred %>% left_join(corrected_symbols_cmg,by = c("GeneName" = "ENSEMBL")) %>% dplyr::rename( gene = SYMBOL)

```


```{r}
gnomAD_v4_constraints_gene <- gnomAD_v4_constraints %>% dplyr::group_by(gene) %>% summarise(lof_hc_lc.pLI.mean = mean(lof_hc_lc.pLI),
                                                                                            lof_hc_lc.pLI.max = max(lof_hc_lc.pLI),
                                                                                            lof_hc_lc.pRec.mean = mean(lof_hc_lc.pRec),
                                                                                            lof_hc_lc.pRec.max = max(lof_hc_lc.pRec),
                                                                                            lof_hc_lc.pNull.mean = max(lof_hc_lc.pNull),
                                                                                            lof_hc_lc.pNull.max = max(lof_hc_lc.pNull),)
Watershed_10k_medz_pred_symbol_gnomad_v4_constraints <- Watershed_10k_medz_pred_symbol %>% left_join(gnomAD_v4_constraints_gene, by = 'gene')
Watershed_10k_cmg_medz_pred_symbol_gnomad_v4_constraints <- Watershed_10k_cmg_medz_pred_symbol %>% left_join(gnomAD_v4_constraints_gene, by = 'gene')

```
```{r}
Watershed_10k_cmg_medz_pred_symbol_gnomad_v4_constraints
```

# mean pLI per gene across transcripts against watershed, gam. 
```{r}
Watershed_10k_medz_pred_symbol_gnomad_v4_constraints %>% 
  filter(!is.na(lof_hc_lc.pNull.mean)) %>% 
  mutate(Watershed_bin = cut(Watershed_posterior_outlier_signal_1, breaks=c(0,0.25,0.5,0.75,1), labels=c("0-0.25","0.25-0.5","0.5-0.75","0.75-1")),
         GAM_bin = cut(GAM_posterior_outlier_signal_1, breaks=c(0,0.25,0.5,0.75,1), labels=c("0-0.25","0.25-0.5","0.5-0.75","0.75-1"))) %>% 
  ggplot(aes(x = Watershed_bin, 
             y = lof_hc_lc.pNull.mean)) + geom_boxplot() #+ scale_color_gradient(low = "blue", high = "red")
```
```{r}
Watershed_10k_cmg_medz_pred_symbol_gnomad_v4_constraints %>% 
  filter(!is.na(lof_hc_lc.pLI.mean)) %>% 
  mutate(Watershed_bin = cut(Watershed_posterior_outlier_signal_1, breaks=c(0,0.25,0.5,0.75,1), labels=c("0-0.25","0.25-0.5","0.5-0.75","0.75-1")),
         GAM_bin = cut(GAM_posterior_outlier_signal_1, breaks=c(0,0.25,0.5,0.75,1), labels=c("0-0.25","0.25-0.5","0.5-0.75","0.75-1"))) %>% 
  ggplot(aes(x = Watershed_bin, 
             y = lof_hc_lc.pLI.mean)) + geom_boxplot() #+ scale_color_gradient(low = "blue", high = "red")
```
```{r}
test <- Watershed_10k_cmg_medz_pred_symbol_gnomad_v4_constraints %>% 
  filter(!is.na(lof_hc_lc.pLI.mean)) %>% 
  mutate(Watershed_bin = cut(Watershed_posterior_outlier_signal_1, breaks=c(0,0.25,0.5,0.75,1), labels=c("0-0.25","0.25-0.5","0.5-0.75","0.75-1")),
         GAM_bin = cut(GAM_posterior_outlier_signal_1, breaks=c(0,0.25,0.5,0.75,1), labels=c("0-0.25","0.25-0.5","0.5-0.75","0.75-1")))
table(test$Watershed_bin)
```

```{r}
Watershed_10k_medz_pred_symbol_gnomad_v4_constraints %>% 
  filter(!is.na(lof_hc_lc.pNull.mean)) %>% 
  mutate(Watershed_bin = cut(Watershed_posterior_outlier_signal_1, breaks=c(0,0.25,0.5,0.75,1), labels=c("0-0.25","0.25-0.5","0.5-0.75","0.75-1")),
         GAM_bin = cut(GAM_posterior_outlier_signal_1, breaks=c(0,0.25,0.5,0.75,1), labels=c("0-0.25","0.25-0.5","0.5-0.75","0.75-1"))) %>% 
  ggplot(aes(x = GAM_bin, 
             y = lof_hc_lc.pNull.mean)) + geom_boxplot() #+ scale_color_gradient(low = "blue", high = "red")
```

```{r}
Watershed_10k_cmg_medz_pred_symbol_gnomad_v4_constraints %>% filter((Watershed_posterior_outlier_signal_1 > 0.5) & (SubjectID=='B14-07')) %>% arrange(desc(Watershed_posterior_outlier_signal_1))
```
```{r}
annotations_cmg <- read.delim('/oak/stanford/groups/smontgom/tannerj/UDN_LRS_watershed/Watershed-SV/MuscularDystrophy-protein-lincRNA-100k-all-position-aware-3.0/MedZ_model_prediction_z3_testing_data.tsv',sep = '\t')
```
```{r}
debug <- annotations_cmg %>% dplyr::select(SVTYPE_DUP_gene_body,
                                  SVTYPE_DUP_tss_flank,
                                  SVTYPE_DUP_tes_flank,
                                  #SVTYPE_DUP_CNV_gene_body,
                                  #SVTYPE_DUP_CNV_tss_flank,
                                  #SVTYPE_DUP_CNV_tes_flank,
                                  SVTYPE_DEL_gene_body,
                                  SVTYPE_DEL_tss_flank,
                                  SVTYPE_DEL_tes_flank,
                                  #SVTYPE_DEL_CNV_gene_body,
                                  #SVTYPE_DEL_CNV_tss_flank,
                                  #SVTYPE_DEL_CNV_tes_flank,
                                  SVTYPE_INS_gene_body,
                                  SVTYPE_INS_tss_flank,
                                  SVTYPE_INS_tes_flank,
                                  SVTYPE_INV_gene_body,
                                  SVTYPE_INV_tss_flank,
                                  SVTYPE_INV_tes_flank) %>% rowSums(na.rm = TRUE)
sum(debug==0)
```

```{r}
annotations_cmg<- annotations_cmg %>% separate(GeneName, into = c("GeneName","SV"), sep = ":", extra = "merge")
```
```{r}
annotations_cmg
```
```{r}
disease_genes <- read_delim(
  'OT-EFO_0002970-associated-targets-1_14_2024-v23_12.tsv',delim = '\t'
)
```
```{r}
disease_genes_confident <- disease_genes %>% filter(globalScore>0.5)
Watershed_10k_cmg_medz_pred_symbol %>% left_join(disease_genes_confident,by = c('gene'='symbol')) %>% mutate(muscle_disease_gene = ifelse(is.na(globalScore),F,T)) %>% filter(muscle_disease_gene==T & Watershed_posterior_outlier_signal_1>0.5)
```
```{r}
Watershed_10k_cmg_medz_pred_symbol_w_annotations %>% filter(GeneName == 'ENSG00000198947' & SubjectID=='251DW_SD')
```


```{r}
Watershed_10k_cmg_medz_pred_symbol_w_annotations <- Watershed_10k_cmg_medz_pred_symbol %>% 
  left_join(disease_genes_confident,by = c('gene'='symbol')) %>% 
  mutate(muscle_disease_gene = ifelse(is.na(globalScore),F,T)) %>% 
  left_join(annotations_cmg,by = c('SubjectID','GeneName','SV')) 
#%>% ggplot(aes(x=-log10(abs(TE_pvalues_+0.0000000001)),y=Watershed_posterior_outlier_signal_1,color = factor(exon_variant))) + geom_point() + ylab('Watershed Posterior') + xlab('-log10(outlier p-value + 1e-10)') +
#  labs(color = 'is exon variant')
```
```{r}
Watershed_10k_cmg_medz_pred_symbol_w_annotations%>% pull('SubjectID') %>% unique()
```

```{r}
posterior_comparison_plot <- Watershed_10k_cmg_medz_pred_symbol %>% 
  left_join(disease_genes_confident,by = c('gene'='symbol')) %>% 
  mutate(muscle_disease_gene = ifelse(is.na(globalScore),F,T)) %>% 
  left_join(annotations_cmg,by = c('SubjectID','GeneName','SV')) %>% 
  mutate(outlier_status = case_when(TE_pvalues_ < -0.0027 ~ "non-outlier",
                                    TE_pvalues_ < 0 ~ "underexpression-outlier",
                                    TE_pvalues_ < 0.0027 ~ "overexpression-outlier",
                                    .default = "non-outlier"),
         variant_category = ifelse(exon_variant==1,'coding','noncoding')) %>% 
  dplyr::select(outlier_status,Watershed_posterior_outlier_signal_1,GAM_posterior_outlier_signal_1,variant_category) %>% 
  pivot_longer(
    cols = c(Watershed_posterior_outlier_signal_1, GAM_posterior_outlier_signal_1), 
    names_to = "Model", 
    values_to = "Posterior"
  ) %>% mutate(Model = ifelse(Model=='Watershed_posterior_outlier_signal_1','Watershed','WGS-only')) %>%
  ggplot(aes(x=outlier_status,y=Posterior,color = factor(Model))) + geom_boxplot() + facet_grid(rows = vars(variant_category), scales="free", space="free_y") + labs(color = 'Model') + xlab('expression outlier status') + ylab('Posterior') + scale_color_manual(values=c("#2C70B0","#FFBC7B"))
posterior_comparison_plot %>% ggsave(file="Figure3/Compare_Watershed_GAM_posterior_CMG.svg", width=7, height=6, dpi=300)

```

```{r}
Watershed_10k_UDN_prediction_file <- '/oak/stanford/groups/smontgom/tannerj/UDN_LRS_watershed/Watershed-SV/UDN-protein-lincRNA-100k-ADRC-all-position-aware-2.0/tissue_specific_annotations_pred_new_impute_SV_Watershed_approximate_100000_number_of_dimensions_2_pthresh_0.0027_trainP_2.56__posterior_probability.txt'
Watershed_10k_UDN_2_tissue_pred <- read_delim(Watershed_10k_UDN_prediction_file,delim = '\t')
Watershed_10k_UDN_2_tissue_pred <- Watershed_10k_UDN_2_tissue_pred %>% separate(sample_names, into = c("SubjectID","GeneName","SV"), sep = ":", extra = "merge")

Watershed_10k_UDN_Medz_blood_prediction_file <- '/oak/stanford/groups/smontgom/tannerj/UDN_LRS_watershed/Watershed-SV/UDN-protein-lincRNA-100k-ADRC-all-position-aware-2.0/MedZ_prediction_z3_Blood_SV_RIVER_100000_number_of_dimensions_1_pthresh_0.0027_trainP_2.56__posterior_probability.txt'
Watershed_10k_UDN_Medz_blood_pred <- read_delim(Watershed_10k_UDN_Medz_blood_prediction_file,delim = '\t')
Watershed_10k_UDN_Medz_blood_pred <- Watershed_10k_UDN_Medz_blood_pred %>% separate(sample_names, into = c("SubjectID","GeneName","SV"), sep = ":", extra = "merge")

Watershed_10k_UDN_Medz_fibroblast_prediction_file <- '/oak/stanford/groups/smontgom/tannerj/UDN_LRS_watershed/Watershed-SV/UDN-protein-lincRNA-100k-ADRC-all-position-aware-2.0/MedZ_prediction_z3_Fibroblast_SV_RIVER_100000_number_of_dimensions_1_pthresh_0.0027_trainP_2.56__posterior_probability.txt'
Watershed_10k_UDN_Medz_fibroblast_pred <- read_delim(Watershed_10k_UDN_Medz_fibroblast_prediction_file,delim = '\t')
Watershed_10k_UDN_Medz_fibroblast_pred <- Watershed_10k_UDN_Medz_fibroblast_pred %>% separate(sample_names, into = c("SubjectID","GeneName","SV"), sep = ":", extra = "merge")
```
```{r}
annotations_udn <- read.delim('/oak/stanford/groups/smontgom/tannerj/UDN_LRS_watershed/Watershed-SV/UDN-protein-lincRNA-100k-ADRC-all-position-aware-2.0/tissue_specific_annotations_pred_no_vep_no_muscle_newimpute_testing_data.tsv',sep = '\t')
annotations_udn_medz_blood <- read.delim('/oak/stanford/groups/smontgom/tannerj/UDN_LRS_watershed/Watershed-SV/UDN-protein-lincRNA-100k-ADRC-all-position-aware-2.0/MedZ_prediction_z3_Blood_testing_data.tsv',sep = '\t')
annotations_udn_medz_fibroblast <- read.delim('/oak/stanford/groups/smontgom/tannerj/UDN_LRS_watershed/Watershed-SV/UDN-protein-lincRNA-100k-ADRC-all-position-aware-2.0/MedZ_prediction_z3_Fibroblast_testing_data.tsv',sep = '\t')
```


```{r}
annotations_udn <- annotations_udn %>% separate(GeneName, into = c("GeneName","SV"), sep = ":", extra = "merge")
annotations_udn_medz_blood <- annotations_udn_medz_blood %>% separate(GeneName, into = c("GeneName","SV"), sep = ":", extra = "merge")
annotations_udn_medz_fibroblast <- annotations_udn_medz_fibroblast %>% separate(GeneName, into = c("GeneName","SV"), sep = ":", extra = "merge")

```

```{r}
posterior_comparison_plot_blood_df <- Watershed_10k_UDN_2_tissue_pred %>% 
  left_join(annotations_udn,by = c('SubjectID','GeneName','SV')) %>%
  filter(!is.na(TE_pvalues_Whole_Blood)) %>%
  mutate(outlier_status = case_when(TE_pvalues_Whole_Blood  < -0.0027 ~ "non-outlier",
                                    TE_pvalues_Whole_Blood < 0 ~ "underexpression",
                                    TE_pvalues_Whole_Blood < 0.0027 ~ "overexpression",
                                    .default = "non-outlier"),
         variant_category = ifelse(exon_variant==1,'coding','noncoding')) %>% 
  dplyr::select(outlier_status,Watershed_posterior_outlier_signal_1,GAM_posterior_outlier_signal_1,variant_category) %>% 
  pivot_longer(
    cols = c(Watershed_posterior_outlier_signal_1, GAM_posterior_outlier_signal_1), 
    names_to = "Model", 
    values_to = "Posterior"
  ) %>% mutate(Model = ifelse(Model=='Watershed_posterior_outlier_signal_1','Watershed','WGS-only')) 

posterior_comparison_plot_blood <- posterior_comparison_plot_blood_df %>%
  ggplot(aes(x=outlier_status,y=Posterior,color = factor(Model))) + geom_boxplot() + facet_grid(rows = vars(variant_category), scales="free", space="free_y") + labs(color = 'Model') + xlab('expression outlier status') + ylab('Posterior') + scale_color_manual(values=c("#2C70B0","#FFBC7B"))
posterior_comparison_plot_blood
```


```{r}
posterior_comparison_plot_blood_medz_df <- Watershed_10k_UDN_Medz_blood_pred %>% 
  left_join(annotations_udn_medz_blood,by = c('SubjectID','GeneName','SV')) %>%
  filter(!is.na(TE_pvalues_)) %>%
  mutate(outlier_status = case_when(TE_pvalues_  < -0.0027 ~ "non-outlier",
                                    TE_pvalues_ < 0 ~ "underexpression",
                                    TE_pvalues_ < 0.0027 ~ "overexpression",
                                    .default = "non-outlier"),
         variant_category = ifelse(exon_variant==1,'coding','noncoding')) %>% 
  dplyr::select(outlier_status,Watershed_posterior_outlier_signal_1,GAM_posterior_outlier_signal_1,variant_category) %>% 
  pivot_longer(
    cols = c(Watershed_posterior_outlier_signal_1, GAM_posterior_outlier_signal_1), 
    names_to = "Model", 
    values_to = "Posterior"
  ) %>% mutate(Model = ifelse(Model=='Watershed_posterior_outlier_signal_1','Watershed','WGS-only')) 
posterior_comparison_plot_blood_medz <- posterior_comparison_plot_blood_medz_df %>%
  ggplot(aes(x=outlier_status,y=Posterior,color = factor(Model))) + geom_boxplot() + facet_grid(rows = vars(variant_category), scales="free", space="free_y") + labs(color = 'Model') + xlab('expression outlier status') + ylab('Posterior') + scale_color_manual(values=c("#2C70B0","#FFBC7B"))
posterior_comparison_plot_blood_medz
```

```{r}
posterior_comparison_plot_fibroblast_df <- Watershed_10k_UDN_2_tissue_pred %>% 
  left_join(annotations_udn,by = c('SubjectID','GeneName','SV')) %>% 
  filter(!is.na(TE_pvalues_Fibroblast)) %>%
  mutate(outlier_status = case_when(TE_pvalues_Fibroblast < -0.0027 ~ "non-outlier",
                                    TE_pvalues_Fibroblast < 0 ~ "underexpression",
                                    TE_pvalues_Fibroblast < 0.0027 ~ "overexpression",
                                    .default = "non-outlier"),
         variant_category = ifelse(exon_variant==1,'coding','noncoding')) %>% 
  dplyr::select(outlier_status,Watershed_posterior_outlier_signal_2,GAM_posterior_outlier_signal_2,variant_category) %>% 
  pivot_longer(
    cols = c(Watershed_posterior_outlier_signal_2, GAM_posterior_outlier_signal_2), 
    names_to = "Model", 
    values_to = "Posterior"
  ) %>% mutate(Model = ifelse(Model=='Watershed_posterior_outlier_signal_2','Watershed','WGS-only')) 
posterior_comparison_plot_fibroblast <- posterior_comparison_plot_fibroblast_df %>%
  ggplot(aes(x=outlier_status,y=Posterior,color = factor(Model))) + geom_boxplot() + facet_grid(rows = vars(variant_category), scales="free", space="free_y") + labs(color = 'Model') + xlab('expression outlier status') + ylab('Posterior') + scale_color_manual(values=c("#2C70B0","#FFBC7B"))
posterior_comparison_plot_fibroblast
```
```{r}
posterior_comparison_plot_fibroblast_medz_df <- Watershed_10k_UDN_Medz_fibroblast_pred %>% 
  left_join(annotations_udn_medz_fibroblast,by = c('SubjectID','GeneName','SV')) %>% 
  filter(!is.na(TE_pvalues_)) %>%
  mutate(outlier_status = case_when(TE_pvalues_ < -0.0027 ~ "non-outlier",
                                    TE_pvalues_ < 0 ~ "underexpression",
                                    TE_pvalues_ < 0.0027 ~ "overexpression",
                                    .default = "non-outlier"),
         variant_category = ifelse(exon_variant==1,'coding','noncoding')) %>% 
  dplyr::select(outlier_status,Watershed_posterior_outlier_signal_1,GAM_posterior_outlier_signal_1,variant_category) %>% 
  pivot_longer(
    cols = c(Watershed_posterior_outlier_signal_1, GAM_posterior_outlier_signal_1), 
    names_to = "Model", 
    values_to = "Posterior"
  ) %>% mutate(Model = ifelse(Model=='Watershed_posterior_outlier_signal_1','Watershed','WGS-only')) 
posterior_comparison_plot_fibroblast_medz <- posterior_comparison_plot_fibroblast_medz_df %>%
  ggplot(aes(x=outlier_status,y=Posterior,color = factor(Model))) + geom_boxplot() + facet_grid(rows = vars(variant_category), scales="free", space="free_y") + labs(color = 'Model') + xlab('expression outlier status') + ylab('Posterior') + scale_color_manual(values=c("#2C70B0","#FFBC7B"))
posterior_comparison_plot_fibroblast_medz
```
#figure 3 supplementary: for other tissues and UDN how does this compare?

```{r}
mutliple_model_posterior_comparison_df <- rbind(posterior_comparison_plot_fibroblast_medz_df %>% mutate(training_RNA_signal = 'median z-score'),
posterior_comparison_plot_blood_medz_df %>% mutate(training_RNA_signal = 'median z-score'),
posterior_comparison_plot_fibroblast_df %>% mutate(training_RNA_signal = 'fibroblast z-score'),
posterior_comparison_plot_blood_df %>% mutate(training_RNA_signal = 'blood z-score'))

grid_coding_noncoding <- mutliple_model_posterior_comparison_df %>% 
  ggplot(aes(x=outlier_status,y=Posterior,color = factor(Model))) + 
  geom_boxplot() + 
  facet_grid(rows = vars(training_RNA_signal), 
             cols = vars(variant_category),
             scales = "free", 
             space= "free_y") + 
  labs(color = 'Model') + 
  xlab('expression outlier status') + 
  ylab('Posterior') + 
  scale_color_manual(values=c("#2C70B0","#FFBC7B"))
grid_coding_noncoding %>% ggsave(file="Figure3/Supplementary_6c_noncoding_coding_by_trained_tissue.svg", width=10, height=7, dpi=300)
```
#figure 3 supplementary: length is important but not trend setting 
```{r}
breaks <- c(10,100,1000,10000,Inf)
labels <- c('>10','>100','>1000','>10000')
posterior_comparison_plot_length <- Watershed_10k_cmg_medz_pred_symbol %>% 
  left_join(disease_genes_confident,by = c('gene'='symbol')) %>% 
  mutate(muscle_disease_gene = ifelse(is.na(globalScore),F,T)) %>% 
  left_join(annotations_cmg,by = c('SubjectID','GeneName','SV')) %>% 
  mutate(outlier_status = case_when(TE_pvalues_ < -0.0027 ~ "non-outlier",
                                    TE_pvalues_ < 0 ~ "underexpression-outlier",
                                    TE_pvalues_ < 0.0027 ~ "overexpression-outlier",
                                    .default = "non-outlier"),
         variant_category = ifelse(exon_variant==1,'coding','noncoding'),
         length_exp = exp(length)) %>% 
  dplyr::select(outlier_status,Watershed_posterior_outlier_signal_1,GAM_posterior_outlier_signal_1,variant_category,length_exp) %>% 
  pivot_longer(
    cols = c(Watershed_posterior_outlier_signal_1, GAM_posterior_outlier_signal_1), 
    names_to = "Model", 
    values_to = "Posterior"
  ) %>% mutate(Model = ifelse(Model=='Watershed_posterior_outlier_signal_1','Watershed','WGS-only'),
               length_bin = cut(length_exp, breaks = breaks, labels = labels)) %>%
  ggplot(aes(x=length_bin,y=Posterior,color = factor(Model))) + geom_boxplot() + labs(color = 'Model') + xlab('length_bins') + ylab('Posterior') + scale_color_manual(values=c("#2C70B0","#FFBC7B"))
```

```{r}
posterior_comparison_plot_length
```
```{r}
posterior_comparison_plot_blood_length_df <- Watershed_10k_UDN_2_tissue_pred %>% 
  left_join(annotations_udn,by = c('SubjectID','GeneName','SV')) %>%
  filter(!is.na(TE_pvalues_Whole_Blood)) %>%
  mutate(outlier_status = case_when(TE_pvalues_Whole_Blood  < -0.0027 ~ "non-outlier",
                                    TE_pvalues_Whole_Blood < 0 ~ "underexpression",
                                    TE_pvalues_Whole_Blood < 0.0027 ~ "overexpression",
                                    .default = "non-outlier"),
         variant_category = ifelse(exon_variant==1,'coding','noncoding'),
         length_exp = exp(length)) %>% 
  dplyr::select(outlier_status,Watershed_posterior_outlier_signal_1,GAM_posterior_outlier_signal_1,variant_category,length_exp) %>% 
  pivot_longer(
    cols = c(Watershed_posterior_outlier_signal_1, GAM_posterior_outlier_signal_1), 
    names_to = "Model", 
    values_to = "Posterior"
  ) %>% mutate(Model = ifelse(Model=='Watershed_posterior_outlier_signal_1','Watershed','WGS-only'),
               length_bin = cut(length_exp, breaks = breaks, labels = labels))

posterior_comparison_plot_blood_length <- posterior_comparison_plot_blood_length_df %>% ggplot(aes(x=length_bin,y=Posterior,color = factor(Model))) + geom_boxplot() + labs(color = 'Model') + xlab('length_bins') + ylab('Posterior') + scale_color_manual(values=c("#2C70B0","#FFBC7B"))

posterior_comparison_plot_blood_length
```


```{r}
posterior_comparison_plot_blood_medz_length_df <- Watershed_10k_UDN_Medz_blood_pred %>% 
  left_join(annotations_udn_medz_blood,by = c('SubjectID','GeneName','SV')) %>%
  filter(!is.na(TE_pvalues_)) %>%
  mutate(outlier_status = case_when(TE_pvalues_  < -0.0027 ~ "non-outlier",
                                    TE_pvalues_ < 0 ~ "underexpression",
                                    TE_pvalues_ < 0.0027 ~ "overexpression",
                                    .default = "non-outlier"),
         variant_category = ifelse(exon_variant==1,'coding','noncoding'),
         length_exp = exp(length)) %>% 
  dplyr::select(outlier_status,Watershed_posterior_outlier_signal_1,GAM_posterior_outlier_signal_1,variant_category,length_exp) %>% 
  pivot_longer(
    cols = c(Watershed_posterior_outlier_signal_1, GAM_posterior_outlier_signal_1), 
    names_to = "Model", 
    values_to = "Posterior"
  ) %>% mutate(Model = ifelse(Model=='Watershed_posterior_outlier_signal_1','Watershed','WGS-only'),
               length_bin = cut(length_exp, breaks = breaks, labels = labels))

posterior_comparison_plot_blood_medz_length <- posterior_comparison_plot_blood_medz_length_df %>% ggplot(aes(x=length_bin,y=Posterior,color = factor(Model))) + geom_boxplot() + labs(color = 'Model') + xlab('length_bins') + ylab('Posterior') + scale_color_manual(values=c("#2C70B0","#FFBC7B"))

posterior_comparison_plot_blood_medz_length 

```
```{r}
posterior_comparison_plot_fibroblast_length_df <- Watershed_10k_UDN_2_tissue_pred %>% 
  left_join(annotations_udn,by = c('SubjectID','GeneName','SV')) %>% 
  filter(!is.na(TE_pvalues_Fibroblast)) %>%
  mutate(outlier_status = case_when(TE_pvalues_Fibroblast < -0.0027 ~ "non-outlier",
                                    TE_pvalues_Fibroblast < 0 ~ "underexpression",
                                    TE_pvalues_Fibroblast < 0.0027 ~ "overexpression",
                                    .default = "non-outlier"),
         variant_category = ifelse(exon_variant==1,'coding','noncoding'),
         length_exp = exp(length)) %>% 
  dplyr::select(outlier_status,Watershed_posterior_outlier_signal_2,GAM_posterior_outlier_signal_2,variant_category,length_exp) %>% 
  pivot_longer(
    cols = c(Watershed_posterior_outlier_signal_2, GAM_posterior_outlier_signal_2), 
    names_to = "Model", 
    values_to = "Posterior"
  ) %>% mutate(Model = ifelse(Model=='Watershed_posterior_outlier_signal_2','Watershed','WGS-only'),
               length_bin = cut(length_exp, breaks = breaks, labels = labels)) 
posterior_comparison_plot_fibroblast_length <- posterior_comparison_plot_fibroblast_length_df %>% ggplot(aes(x=length_bin,y=Posterior,color = factor(Model))) + geom_boxplot() + labs(color = 'Model') + xlab('length_bins') + ylab('Posterior') + scale_color_manual(values=c("#2C70B0","#FFBC7B"))# + facet_grid(rows = vars(outlier_status), scales="free", space="free_y")
posterior_comparison_plot_fibroblast_length
```
```{r}
posterior_comparison_plot_fibroblast_medz_length_df <- Watershed_10k_UDN_Medz_fibroblast_pred %>% 
  left_join(annotations_udn_medz_fibroblast,by = c('SubjectID','GeneName','SV')) %>% 
  filter(!is.na(TE_pvalues_)) %>%
  mutate(outlier_status = case_when(TE_pvalues_ < -0.0027 ~ "non-outlier",
                                    TE_pvalues_ < 0 ~ "underexpression",
                                    TE_pvalues_ < 0.0027 ~ "overexpression",
                                    .default = "non-outlier"),
         variant_category = ifelse(exon_variant==1,'coding','noncoding'),
         length_exp = exp(length)) %>% 
  dplyr::select(outlier_status,Watershed_posterior_outlier_signal_1,GAM_posterior_outlier_signal_1,variant_category,length_exp) %>% 
  pivot_longer(
    cols = c(Watershed_posterior_outlier_signal_1, GAM_posterior_outlier_signal_1), 
    names_to = "Model", 
    values_to = "Posterior"
  ) %>% mutate(Model = ifelse(Model=='Watershed_posterior_outlier_signal_1','Watershed','WGS-only'),
               length_bin = cut(length_exp, breaks = breaks, labels = labels)) 
posterior_comparison_plot_fibroblast_medz_length <- posterior_comparison_plot_fibroblast_medz_length_df %>%
  ggplot(aes(x=length_bin,y=Posterior,color = factor(Model))) + geom_boxplot() + labs(color = 'Model') + xlab('length_bins') + ylab('Posterior') + scale_color_manual(values=c("#2C70B0","#FFBC7B"))# + facet_grid(rows = vars(outlier_status), scales="free", space="free_y")
posterior_comparison_plot_fibroblast_medz_length
```

```{r}
mutliple_model_posterior_comparison_length_df <- 
  rbind(
    posterior_comparison_plot_fibroblast_medz_length_df %>% mutate(training_RNA_signal = 'median z-score',
                                                                   testing_RNA_signal = 'fibroblast'),
    posterior_comparison_plot_blood_medz_length_df %>% mutate(training_RNA_signal = 'median z-score',
                                                              testing_RNA_signal = 'blood'),
    posterior_comparison_plot_fibroblast_length_df %>% mutate(training_RNA_signal = 'tissue z-score',
                                                              testing_RNA_signal = 'fibroblast'),
    posterior_comparison_plot_blood_length_df %>% mutate(training_RNA_signal = 'tissue z-score',
                                                         testing_RNA_signal = 'blood'))

grid_coding_noncoding <- mutliple_model_posterior_comparison_length_df %>% 
  ggplot(aes(x=length_bin,y=Posterior,color = factor(Model))) +
  geom_boxplot() + 
  xlab('length_bins') + 
  ylab('Posterior') + 
  labs(color = 'Model') +
  facet_grid(rows = vars(testing_RNA_signal), 
             cols = vars(training_RNA_signal),
             scales = "free", 
             space= "free_y") + 
  scale_color_manual(values=c("#2C70B0","#FFBC7B"))
grid_coding_noncoding %>% ggsave(file="Figure3/Supplementary_6d_length_vs_posterior_supplementary_by_trained_tissue.svg", width=10, height=7, dpi=300)
```

# this part is reviewer's response. 
# part 1 about all sv. 
```{r}
gtex_watershed_annot <- read_delim('/oak/stanford/groups/smontgom/tannerj/UDN_LRS_watershed/Watershed-SV/new-protein-lincRNA-100k-rare-position-aware-2.0/combined_annotation_pre_merge_noimpute.medZ.csv',
                                delim = '\t')
gtex_caddsv_score <- read_delim('/oak/stanford/groups/smontgom/tannerj/UDN_LRS_watershed/Watershed-SV/input/gtex_score.txt',
                                delim = '\t')
sv_gtex <- read_delim('/oak/stanford/groups/smontgom/tannerj/UDN_LRS_watershed/Watershed-SV/new-protein-lincRNA-100k-rare-position-aware-2.0/intermediates/pipeline_input.bed', 
                      delim = '\t', col_names = c('chrom','start','end','SV','SVTYPE'))

sv_gtex_reformatted <- sv_gtex %>% mutate(chr = str_remove(chrom, "chr"), 
                                          type = ifelse(str_detect(SVTYPE,'_CNV'), 
                                                        str_remove(SVTYPE,'_CNV'),SVTYPE)) %>% dplyr::select(-c('chrom','SVTYPE')) %>% dplyr::select(c('chr','start','end','SV','type'))

gtex_caddsv_score_reformatted <- gtex_caddsv_score %>% 
  left_join(sv_gtex_reformatted, by = c('chr','start','end','type')) %>% 
  dplyr::select(
    1:5,  # Select the first 6 columns
    last_col(),  # Select the last column
    everything()  # Select all remaining columns
  )
```

```{r}
gtex_SV_gene_mapping <- gtex_watershed_annot %>% dplyr::select(SubjectID, GeneName) %>% separate(GeneName, into = c('GeneName','SV'), sep = ':', remove = TRUE)

gtex_caddsv_score_reformatted_reduced <- gtex_caddsv_score_reformatted %>% dplyr::select(c('SV','CADDSV-score'))

gtex_caddsv_score_gene_level <- gtex_SV_gene_mapping %>% 
  left_join(gtex_caddsv_score_reformatted_reduced,by='SV') %>% 
  drop_na(`CADDSV-score`) %>% 
  group_by(SubjectID,GeneName) %>% 
  summarise(CADDSV_gene_score = max(`CADDSV-score`, na.rm = TRUE))

gtex_caddsv_score_gene_level %>% write_delim('/oak/stanford/groups/smontgom/tannerj/UDN_LRS_watershed/Watershed-SV/input/CADDSV_GTEx_prediction_max_cisSVs.tsv', delim = '\t')
```

# noncoding parallel.
```{r}
watershed_ABC_100k_obj_Cells_Cultured_Fibroblasts_file  <- '/oak/stanford/groups/smontgom/tannerj/UDN_LRS_watershed/Watershed-SV/new-protein-lincRNA-100k-rare-position-aware-2.0/eval_watarshed_ABC_regional_Cells_Cultured_fibroblast_SV_RIVER_100000_number_of_dimensions_1_pthresh_0.000063_pval_frac_0.176_separate_trainP_2.56__evaluation_object.rds'
data_Cells_Cultured_Fibroblasts_file <- '/oak/stanford/groups/smontgom/tannerj/UDN_LRS_watershed/Watershed-SV/new-protein-lincRNA-100k-rare-position-aware-2.0/eval_watarshed_ABC_regional.Cells_Cultured_fibroblasts.z4.tsv'
data_Cells_Cultured_Fibroblasts_sv_file <- '/oak/stanford/groups/smontgom/tannerj/UDN_LRS_watershed/Watershed-SV/new-protein-lincRNA-100k-rare-position-aware-2.0/combined_annotation_pre_merge_noimpute.Cells_Cultured_fibroblasts.z4.jointcalled.csv'

watershed_ABC_100k_obj_Muscle_Skeletal_file <- '/oak/stanford/groups/smontgom/tannerj/UDN_LRS_watershed/Watershed-SV/new-protein-lincRNA-100k-rare-position-aware-2.0/eval_watarshed_ABC_regional_Muscle_Skeletal_SV_RIVER_100000_number_of_dimensions_1_pthresh_0.000063_pval_frac_0.151_separate_trainP_2.56__evaluation_object.rds'
data_Muscle_Skeletal_file <- '/oak/stanford/groups/smontgom/tannerj/UDN_LRS_watershed/Watershed-SV/new-protein-lincRNA-100k-rare-position-aware-2.0/eval_watarshed_ABC_regional.Muscle_Skeletal.z4.tsv'
data_Muscle_Skeletal_sv_file <- '/oak/stanford/groups/smontgom/tannerj/UDN_LRS_watershed/Watershed-SV/new-protein-lincRNA-100k-rare-position-aware-2.0/combined_annotation_pre_merge_noimpute.Muscle_Skeletal.z4.jointcalled.csv'

watershed_ABC_100k_obj_Whole_Blood_file <- '/oak/stanford/groups/smontgom/tannerj/UDN_LRS_watershed/Watershed-SV/new-protein-lincRNA-100k-rare-position-aware-2.0/eval_watarshed_ABC_regional_Whole_Blood_SV_RIVER_100000_number_of_dimensions_1_pthresh_0.000063_pval_frac_0.138_separate_trainP_2.56__evaluation_object.rds'
data_Whole_Blood_file <- '/oak/stanford/groups/smontgom/tannerj/UDN_LRS_watershed/Watershed-SV/new-protein-lincRNA-100k-rare-position-aware-2.0/eval_watarshed_ABC_regional.Whole_Blood.z4.tsv'
data_Whole_Blood_sv_file <- '/oak/stanford/groups/smontgom/tannerj/UDN_LRS_watershed/Watershed-SV/new-protein-lincRNA-100k-rare-position-aware-2.0/combined_annotation_pre_merge_noimpute.Whole_Blood.z4.jointcalled.csv'

watershed_ABC_100k_obj_Median_file  <- '/oak/stanford/groups/smontgom/tannerj/UDN_LRS_watershed/Watershed-SV/new-protein-lincRNA-100k-rare-position-aware-2.0/eval_watarshed_ABC_regional_MedZ_SV_RIVER_100000_number_of_dimensions_1_pthresh_0.0027_pval_frac_0.16_separate_trainP_2.56__evaluation_object.rds'
data_Median_file <- '/oak/stanford/groups/smontgom/tannerj/UDN_LRS_watershed/Watershed-SV/new-protein-lincRNA-100k-rare-position-aware-2.0/eval_watarshed_ABC_regional.medZ.z3.tsv'
data_Median_sv_file <- '/oak/stanford/groups/smontgom/tannerj/UDN_LRS_watershed/Watershed-SV/new-protein-lincRNA-100k-rare-position-aware-2.0/combined_annotation_pre_merge_noimpute.medZ.csv'
```

# fibroblast
```{r}
gtex_watershed_annot_fibroblast <- read_delim(data_Cells_Cultured_Fibroblasts_sv_file,
                                delim = '\t')
gtex_caddsv_score <- read_delim('/oak/stanford/groups/smontgom/tannerj/UDN_LRS_watershed/Watershed-SV/input/gtex_score.txt',
                                delim = '\t')
sv_gtex <- read_delim('/oak/stanford/groups/smontgom/tannerj/UDN_LRS_watershed/Watershed-SV/new-protein-lincRNA-100k-rare-position-aware-2.0/intermediates/pipeline_input.bed', 
                      delim = '\t', col_names = c('chrom','start','end','SV','SVTYPE'))

sv_gtex_reformatted <- sv_gtex %>% mutate(chr = str_remove(chrom, "chr"), 
                                          type = ifelse(str_detect(SVTYPE,'_CNV'), 
                                                        str_remove(SVTYPE,'_CNV'),SVTYPE)) %>% dplyr::select(-c('chrom','SVTYPE')) %>% dplyr::select(c('chr','start','end','SV','type'))

gtex_caddsv_score_reformatted <- gtex_caddsv_score %>% 
  left_join(sv_gtex_reformatted, by = c('chr','start','end','type')) %>% 
  dplyr::select(
    1:5,  # Select the first 6 columns
    last_col(),  # Select the last column
    everything()  # Select all remaining columns
  )

gtex_SV_gene_mapping_fibroblast <- gtex_watershed_annot_fibroblast %>% dplyr::select(SubjectID, GeneName) %>% separate(GeneName, into = c('GeneName','SV'), sep = ':', remove = TRUE)

gtex_caddsv_score_reformatted_reduced <- gtex_caddsv_score_reformatted %>% dplyr::select(c('SV','CADDSV-score'))

gtex_caddsv_score_gene_level_fibroblast <- gtex_SV_gene_mapping_fibroblast %>% 
  left_join(gtex_caddsv_score_reformatted_reduced,by='SV') %>% 
  drop_na(`CADDSV-score`) %>% 
  group_by(SubjectID,GeneName) %>% 
  summarise(CADDSV_gene_score = max(`CADDSV-score`, na.rm = TRUE))

gtex_caddsv_score_gene_level_fibroblast %>% write_delim('/oak/stanford/groups/smontgom/tannerj/UDN_LRS_watershed/Watershed-SV/input/CADDSV_GTEx_prediction_max_cisSVs.Fibroblast.tsv', delim = '\t')
```
# muscle
```{r}
gtex_watershed_annot_muscle <- read_delim(data_Muscle_Skeletal_sv_file,
                                delim = '\t')
gtex_caddsv_score <- read_delim('/oak/stanford/groups/smontgom/tannerj/UDN_LRS_watershed/Watershed-SV/input/gtex_score.txt',
                                delim = '\t')
sv_gtex <- read_delim('/oak/stanford/groups/smontgom/tannerj/UDN_LRS_watershed/Watershed-SV/new-protein-lincRNA-100k-rare-position-aware-2.0/intermediates/pipeline_input.bed', 
                      delim = '\t', col_names = c('chrom','start','end','SV','SVTYPE'))

sv_gtex_reformatted <- sv_gtex %>% mutate(chr = str_remove(chrom, "chr"), 
                                          type = ifelse(str_detect(SVTYPE,'_CNV'), 
                                                        str_remove(SVTYPE,'_CNV'),SVTYPE)) %>% dplyr::select(-c('chrom','SVTYPE')) %>% dplyr::select(c('chr','start','end','SV','type'))

gtex_caddsv_score_reformatted <- gtex_caddsv_score %>% 
  left_join(sv_gtex_reformatted, by = c('chr','start','end','type')) %>% 
  dplyr::select(
    1:5,  # Select the first 6 columns
    last_col(),  # Select the last column
    everything()  # Select all remaining columns
  )

gtex_SV_gene_mapping_muscle <- gtex_watershed_annot_muscle %>% dplyr::select(SubjectID, GeneName) %>% separate(GeneName, into = c('GeneName','SV'), sep = ':', remove = TRUE)

gtex_caddsv_score_reformatted_reduced <- gtex_caddsv_score_reformatted %>% dplyr::select(c('SV','CADDSV-score'))

gtex_caddsv_score_gene_level_muscle <- gtex_SV_gene_mapping_muscle %>% 
  left_join(gtex_caddsv_score_reformatted_reduced,by='SV') %>% 
  drop_na(`CADDSV-score`) %>% 
  group_by(SubjectID,GeneName) %>% 
  summarise(CADDSV_gene_score = max(`CADDSV-score`, na.rm = TRUE))

gtex_caddsv_score_gene_level_muscle %>% write_delim('/oak/stanford/groups/smontgom/tannerj/UDN_LRS_watershed/Watershed-SV/input/CADDSV_GTEx_prediction_max_cisSVs.Muscle.tsv', delim = '\t')
```

# blood
```{r}
gtex_watershed_annot_blood <- read_delim(data_Whole_Blood_sv_file,
                                delim = '\t')
gtex_caddsv_score <- read_delim('/oak/stanford/groups/smontgom/tannerj/UDN_LRS_watershed/Watershed-SV/input/gtex_score.txt',
                                delim = '\t')
sv_gtex <- read_delim('/oak/stanford/groups/smontgom/tannerj/UDN_LRS_watershed/Watershed-SV/new-protein-lincRNA-100k-rare-position-aware-2.0/intermediates/pipeline_input.bed', 
                      delim = '\t', col_names = c('chrom','start','end','SV','SVTYPE'))

sv_gtex_reformatted <- sv_gtex %>% mutate(chr = str_remove(chrom, "chr"), 
                                          type = ifelse(str_detect(SVTYPE,'_CNV'), 
                                                        str_remove(SVTYPE,'_CNV'),SVTYPE)) %>% dplyr::select(-c('chrom','SVTYPE')) %>% dplyr::select(c('chr','start','end','SV','type'))

gtex_caddsv_score_reformatted <- gtex_caddsv_score %>% 
  left_join(sv_gtex_reformatted, by = c('chr','start','end','type')) %>% 
  dplyr::select(
    1:5,  # Select the first 6 columns
    last_col(),  # Select the last column
    everything()  # Select all remaining columns
  )

gtex_SV_gene_mapping_blood <- gtex_watershed_annot_blood %>% dplyr::select(SubjectID, GeneName) %>% separate(GeneName, into = c('GeneName','SV'), sep = ':', remove = TRUE)

gtex_caddsv_score_reformatted_reduced <- gtex_caddsv_score_reformatted %>% dplyr::select(c('SV','CADDSV-score'))

gtex_caddsv_score_gene_level_blood <- gtex_SV_gene_mapping_blood %>% 
  left_join(gtex_caddsv_score_reformatted_reduced,by='SV') %>% 
  drop_na(`CADDSV-score`) %>% 
  group_by(SubjectID,GeneName) %>% 
  summarise(CADDSV_gene_score = max(`CADDSV-score`, na.rm = TRUE))

gtex_caddsv_score_gene_level_blood %>% write_delim('/oak/stanford/groups/smontgom/tannerj/UDN_LRS_watershed/Watershed-SV/input/CADDSV_GTEx_prediction_max_cisSVs.Blood.tsv', delim = '\t')
```
```{r}

```

