---
title: "Revision_CADDSV_WatershedSV"
output: html_document
date: '2024-08-15'
---

```{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}
remove_na <- function(x) {
	return(x[!is.na(x)])
}
# this version updated so can compute caddsv comparison at the same time. 
compute_auc_across_dimensions_caddsv <- function(number_of_dimensions, watershed_posteriors, gam_posteriors, caddsv_score, real_valued_outliers_test1, binary_outliers_test2) {
	# Initialize objects to store auc information
	auc_object_across_dimensions <- list()
	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)
      
   		# CADDSV evaluation curves
  		caddsv_roc_obj <- roc.curve(scores.class0 = remove_na(caddsv_score[,dimension][test_outlier_status==1 & !is.na(real_valued_outliers_test1[,dimension])]), scores.class1 = remove_na(caddsv_score[,dimension][test_outlier_status==0 & !is.na(real_valued_outliers_test1[,dimension])]), curve = T)
   		caddsv_pr_obj <- pr.curve(scores.class0 = remove_na(caddsv_score[,dimension][test_outlier_status==1 & !is.na(real_valued_outliers_test1[,dimension])]), scores.class1 = remove_na(caddsv_score[,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],
		      CADDSV_sens=caddsv_roc_obj$curve[,2],
          CADDSV_spec=1-caddsv_roc_obj$curve[,1],
          CADDSV_auc=caddsv_roc_obj$auc,
         	CADDSV_pr_auc=caddsv_pr_obj$auc.integral,
         	CADDSV_recall=caddsv_pr_obj$curve[,1],
         	CADDSV_precision=caddsv_pr_obj$curve[,2])

		auc_object_across_dimensions[[dimension]] <- list(evaROC=evaROC)
	}
	return(auc_object_across_dimensions)
}

re_calculate_precision_recall_caddsv <- function(watershed_object, eval_data,tissue_count,fraction_outlier,pvalue_threshold,compute_ci = FALSE) {
  watershed_ABC_100k_obj <- 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
  cadd <- watershed_ABC_100k_obj$n2_posteriors$caddsv_n2_score
  rvot1 <- real_valued_outliers_test1
  bot2 <- binary_outliers_test2_alt
  auc_object_across_dimensions <- compute_auc_across_dimensions_caddsv(1,
                                                                watershed,
                                                                gam, 
                                                                cadd,
                                                                rvot1, 
                                                                bot2)
  return (auc_object_across_dimensions)
  }
}
re_calculate_precision_recall_caddsv_noncoding <- function(watershed_object,eval_data,tissue_count,fraction_outlier,pvalue_threshold,compute_ci = FALSE) {
  watershed_ABC_100k_obj <- 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),])
    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 <- 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),])
    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],
              caddsv_n2_score = watershed_ABC_100k_obj$n2_posteriors$caddsv_n2_score[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_caddsv <- matrix(watershed_ABC_100k_obj$n2_posteriors$caddsv_n2_score[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_caddsv(1,
                                                              matrix(noncoding_watershed),
                                                              matrix(noncoding_gam), 
                                                              matrix(noncoding_caddsv), 
                                                              matrix(noncoding_rvot1), 
                                                              noncoding_bot2)
    return(auc_object_across_dimensions)
  }
}
```
```{r}
# debug block
watershed_load_debug <- load_watershed_data(data_Median_file,1,0.176,0.05)
```

```{r}
feat_all <- watershed_load_debug$feat
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),])

```
```{r}
data.frame(feat_test)
```

```{r}
re_calculate_precision_recall_caddsv_by_distance <- function(watershed_object,eval_data, tissue_count,fraction_outlier,pvalue_threshold, distance_bin = 'intragenic coding', compute_ci = FALSE) {
  watershed_ABC_100k_obj <- 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)),]
    if ( distance_bin == 'intragenic coding or noncoding'){
      is_in_distance_category_1 <- c((feat_test[,"nc_dist2gene" ]==0)[seq(from=1,to=sum((!is.na(N2_pairs))),by=2)],
                                     (feat_test[,"nc_dist2gene" ]==0)[seq(from=2,to=sum((!is.na(N2_pairs))),by=2)])
      
      is_in_distance_category_2 <- c((feat_test[,"nc_dist2gene" ]==0)[seq(from=2,to=sum((!is.na(N2_pairs))),by=2)],
                                     (feat_test[,"nc_dist2gene" ]==0)[seq(from=1,to=sum((!is.na(N2_pairs))),by=2)])

    } else if ( distance_bin == '1b-50kb noncoding' ){
      is_in_distance_category_1 <- c((feat_test[,"coding_fraction_affected" ]==0 & feat_test[,"nc_dist2gene" ] > 0 & feat_test[,"nc_dist2gene" ] <=50000)[seq(from=1,to=sum((!is.na(N2_pairs))),by=2)],
                  (feat_test[,"coding_fraction_affected" ]==0 & feat_test[,"nc_dist2gene" ] > 0 & feat_test[,"nc_dist2gene" ] <=50000)[seq(from=2,to=sum((!is.na(N2_pairs))),by=2)])
      is_in_distance_category_2 <- c((feat_test[,"coding_fraction_affected" ]==0 & feat_test[,"nc_dist2gene" ] > 0 & feat_test[,"nc_dist2gene" ] <=50000)[seq(from=2,to=sum((!is.na(N2_pairs))),by=2)],
                  (feat_test[,"coding_fraction_affected" ]==0 & feat_test[,"nc_dist2gene" ] > 0 & feat_test[,"nc_dist2gene" ] <=50000)[seq(from=1,to=sum((!is.na(N2_pairs))),by=2)])
    } else if ( distance_bin == '50kb-100kb noncoding' ){
      is_in_distance_category_1 <- c((feat_test[,"coding_fraction_affected" ]==0 & feat_test[,"nc_dist2gene" ] > 50000 & feat_test[,"nc_dist2gene" ] <=100000)[seq(from=1,to=sum((!is.na(N2_pairs))),by=2)],
                  (feat_test[,"coding_fraction_affected" ]==0 & feat_test[,"nc_dist2gene" ] > 50000 & feat_test[,"nc_dist2gene" ] <=100000)[seq(from=2,to=sum((!is.na(N2_pairs))),by=2)])
      is_in_distance_category_2 <- c((feat_test[,"coding_fraction_affected" ]==0 & feat_test[,"nc_dist2gene" ] > 50000 & feat_test[,"nc_dist2gene" ] <=100000)[seq(from=2,to=sum((!is.na(N2_pairs))),by=2)],
                  (feat_test[,"coding_fraction_affected" ]==0 & feat_test[,"nc_dist2gene" ] > 50000 & feat_test[,"nc_dist2gene" ] <=100000)[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),]
    if ( distance_bin == 'intragenic coding or noncoding' ){
      is_in_distance_category_1 <- c((feat_test[,"nc_dist2gene" ]==0)[seq(from=1,to=sum((!is.na(N2_pairs))),by=2)],
                  (feat_test[,"nc_dist2gene" ]==0)[seq(from=2,to=sum((!is.na(N2_pairs))),by=2)])
      is_in_distance_category_2 <- c((feat_test[,"nc_dist2gene" ]==0)[seq(from=2,to=sum((!is.na(N2_pairs))),by=2)],
                  (feat_test[,"nc_dist2gene" ]==0)[seq(from=1,to=sum((!is.na(N2_pairs))),by=2)])
    } else if ( distance_bin == '1b-50kb noncoding' ){
      is_in_distance_category_1 <- c((feat_test[,"coding_fraction_affected" ]==0 & feat_test[,"nc_dist2gene" ] > 0 & feat_test[,"nc_dist2gene" ] <=50000)[seq(from=1,to=sum((!is.na(N2_pairs))),by=2)],
                  (feat_test[,"coding_fraction_affected" ]==0 & feat_test[,"nc_dist2gene" ] > 0 & feat_test[,"nc_dist2gene" ] <=50000)[seq(from=2,to=sum((!is.na(N2_pairs))),by=2)])
      is_in_distance_category_2 <- c((feat_test[,"coding_fraction_affected" ]==0 & feat_test[,"nc_dist2gene" ] > 0 & feat_test[,"nc_dist2gene" ] <=50000)[seq(from=2,to=sum((!is.na(N2_pairs))),by=2)],
                  (feat_test[,"coding_fraction_affected" ]==0 & feat_test[,"nc_dist2gene" ] > 0 & feat_test[,"nc_dist2gene" ] <=50000)[seq(from=1,to=sum((!is.na(N2_pairs))),by=2)])
    } else if ( distance_bin == '50kb-100kb noncoding' ){
      is_in_distance_category_1 <- c((feat_test[,"coding_fraction_affected" ]==0 & feat_test[,"nc_dist2gene" ] > 50000 & feat_test[,"nc_dist2gene" ] <=100000)[seq(from=1,to=sum((!is.na(N2_pairs))),by=2)],
                  (feat_test[,"coding_fraction_affected" ]==0 & feat_test[,"nc_dist2gene" ] > 50000 & feat_test[,"nc_dist2gene" ] <=100000)[seq(from=2,to=sum((!is.na(N2_pairs))),by=2)])
      is_in_distance_category_2 <- c((feat_test[,"coding_fraction_affected" ]==0 & feat_test[,"nc_dist2gene" ] > 50000 & feat_test[,"nc_dist2gene" ] <=100000)[seq(from=2,to=sum((!is.na(N2_pairs))),by=2)],
                  (feat_test[,"coding_fraction_affected" ]==0 & feat_test[,"nc_dist2gene" ] > 50000 & feat_test[,"nc_dist2gene" ] <=100000)[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_in_distance_category_1,1],
              gam_n2_posterior = watershed_ABC_100k_obj$n2_posteriors$gam_n2_posterior[is_in_distance_category_1,1],
              caddsv_n2_score = watershed_ABC_100k_obj$n2_posteriors$caddsv_n2_score[is_in_distance_category_1,1],
              rvot1 = real_valued_outliers_test1[is_in_distance_category_1,1],
              bot2 = binary_outliers_test2_alt[is_in_distance_category_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_in_distance_category_1,])
    noncoding_gam <- matrix(watershed_ABC_100k_obj$n2_posteriors$gam_n2_posterior[is_in_distance_category_1,])
    noncoding_caddsv <- matrix(watershed_ABC_100k_obj$n2_posteriors$caddsv_n2_score[is_in_distance_category_1,])
    noncoding_rvot1 <- matrix(real_valued_outliers_test1[is_in_distance_category_1,])
    noncoding_bot2 <- matrix(binary_outliers_test2_alt[is_in_distance_category_2,])
    auc_object_across_dimensions <- compute_auc_across_dimensions_caddsv(1,
                                                              noncoding_watershed,
                                                              noncoding_gam, 
                                                              noncoding_caddsv, 
                                                              noncoding_rvot1, 
                                                              noncoding_bot2)
    return(auc_object_across_dimensions)
  }
}

make_panel_4d_caddsv <- function(df) {
    # Set ordering of model types
	df$prediction_type=factor(df$prediction_type, levels=c("Watershed", 
	                                                       "WGS-only",
	                                                       "CADDSV"))

    # 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","#D81B60")) +
                scale_linetype_manual(values=c("solid", "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","#D81B60")) +
                scale_linetype_manual(values=c("solid", "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","#D81B60")) +
                scale_linetype_manual(values=c("solid", "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","#D81B60")) +
                scale_linetype_manual(values=c("solid", "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)
}

```

# 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'
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'
data_Median_10k_sv_file <- '/oak/stanford/groups/smontgom/tannerj/UDN_LRS_watershed/Watershed-SV/new-protein-lincRNA-10k-final-version/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')
```

# median
```{r}
gtex_watershed_annot_med <- read_delim(data_Median_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_med <- gtex_watershed_annot_med %>% 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_med <- gtex_SV_gene_mapping_med %>% 
  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_med %>% write_delim('/oak/stanford/groups/smontgom/tannerj/UDN_LRS_watershed/Watershed-SV/input/CADDSV_GTEx_prediction_max_cisSVs.MedZ.tsv', delim = '\t')
```

# median
```{r}
gtex_watershed_annot_med_10k <- read_delim(data_Median_10k_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-10k-final-version/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_med_10k <- gtex_watershed_annot_med_10k %>% 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_med_10k <- gtex_SV_gene_mapping_med_10k %>% 
  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_med_10k %>% write_delim('/oak/stanford/groups/smontgom/tannerj/UDN_LRS_watershed/Watershed-SV/input/CADDSV_GTEx_prediction_max_cisSVs.10k.MedZ.tsv', delim = '\t')
```


# all data filenames
```{r}
gtex_caddsv_score_gene_level <- read_delim('/oak/stanford/groups/smontgom/tannerj/UDN_LRS_watershed/Watershed-SV/input/CADDSV_GTEx_prediction_max_cisSVs.tsv', delim = '\t')
gtex_caddsv_score_gene_level_fibroblast <- read_delim('/oak/stanford/groups/smontgom/tannerj/UDN_LRS_watershed/Watershed-SV/input/CADDSV_GTEx_prediction_max_cisSVs.Fibroblast.tsv', delim = '\t')
gtex_caddsv_score_gene_level_muscle <- read_delim('/oak/stanford/groups/smontgom/tannerj/UDN_LRS_watershed/Watershed-SV/input/CADDSV_GTEx_prediction_max_cisSVs.Muscle.tsv', delim = '\t')
gtex_caddsv_score_gene_level_blood <- read_delim('/oak/stanford/groups/smontgom/tannerj/UDN_LRS_watershed/Watershed-SV/input/CADDSV_GTEx_prediction_max_cisSVs.Blood.tsv', delim = '\t')


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'
```


```{r}
testobj <- readRDS(watershed_ABC_100k_obj_Median_file)
test_features <- testobj$feat_test %>% as.data.frame() %>% rownames_to_column(var = "ID")
cadd_score_pre_merge <- gtex_caddsv_score_gene_level %>%  mutate(ID = paste(SubjectID, GeneName, sep = ":")) %>% dplyr::select(ID,CADDSV_gene_score)
testobj$n2_posteriors$caddsv_n2_score <- test_features %>% dplyr::select(ID) %>% left_join(cadd_score_pre_merge, by = "ID") %>% dplyr::pull(CADDSV_gene_score) %>% matrix(ncol = 1)

testobj_10k <- readRDS(watershed_ABC_10k_obj_Median_file)
test_features_10k <- testobj_10k$feat_test %>% as.data.frame() %>% rownames_to_column(var = "ID")
cadd_score_pre_merge <- gtex_caddsv_score_gene_level %>%  mutate(ID = paste(SubjectID, GeneName, sep = ":")) %>% dplyr::select(ID,CADDSV_gene_score)
testobj_10k$n2_posteriors$caddsv_n2_score <- test_features_10k %>% dplyr::select(ID) %>% left_join(cadd_score_pre_merge, by = "ID") %>% dplyr::pull(CADDSV_gene_score) %>% matrix(ncol = 1)
```


# test if this works
```{r}
median_100k_prec_rec <- re_calculate_precision_recall_caddsv(testobj,data_Median_file,1,0.176,0.05,compute_ci = F)
median_10k_prec_rec <- re_calculate_precision_recall_caddsv(testobj_10k,data_Median_10k_file,1,0.31,0.05,compute_ci = F)

```
```{r}
tissue_dim <- 1
medz_watershed_100kb <- data_frame(recall = median_100k_prec_rec[[tissue_dim]]$evaROC$watershed_recall,
           precision = median_100k_prec_rec[[tissue_dim]]$evaROC$watershed_precision,
           prediction_type = 'Watershed-100kb',outlier_type = 'Median') %>% arrange(desc(row_number()))

medz_gam_100kb <- data_frame(recall = median_100k_prec_rec[[tissue_dim]]$evaROC$GAM_recall,
           precision = median_100k_prec_rec[[tissue_dim]]$evaROC$GAM_precision,
           prediction_type = 'WGS-only-100kb',outlier_type = 'Median') %>% arrange(desc(row_number()))

medz_caddsv_100kb <- data_frame(recall = median_100k_prec_rec[[tissue_dim]]$evaROC$CADDSV_recall,
           precision = median_100k_prec_rec[[tissue_dim]]$evaROC$CADDSV_precision,
           prediction_type = 'CADDSV-100kb',outlier_type = 'Median') %>% arrange(desc(row_number()))

medz_watershed_10kb <- data_frame(recall = median_10k_prec_rec[[tissue_dim]]$evaROC$watershed_recall,
           precision = median_10k_prec_rec[[tissue_dim]]$evaROC$watershed_precision,
           prediction_type = 'Watershed-10kb',outlier_type = 'Median') %>% arrange(desc(row_number()))

medz_gam_10kb <- data_frame(recall = median_10k_prec_rec[[tissue_dim]]$evaROC$GAM_recall,
           precision = median_10k_prec_rec[[tissue_dim]]$evaROC$GAM_precision,
           prediction_type = 'WGS-only-10kb',outlier_type = 'Median') %>% arrange(desc(row_number()))

medz_caddsv_10kb <- data_frame(recall = median_10k_prec_rec[[tissue_dim]]$evaROC$CADDSV_recall,
           precision = median_10k_prec_rec[[tissue_dim]]$evaROC$CADDSV_precision,
           prediction_type = 'CADDSV-10kb',outlier_type = 'Median') %>% arrange(desc(row_number()))

total_medz_pr_data <- rbind(medz_watershed_100kb,medz_gam_100kb,medz_caddsv_100kb,
                            medz_watershed_10kb,medz_gam_10kb,medz_caddsv_10kb)
```

```{r}
make_panel_medz_w_cadd <- 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",
	                                                       "CADDSV-10kb",
	                                                       "CADDSV-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","#D81B60","#D81B60")) +
                scale_linetype_manual(values=c("solid","dashed","solid","dashed","solid","dashed")) +
                gtex_v8_figure_theme() + 
                ggtitle(outlier_name,)
  	return (plotter_ase)
}
```
```{r}
p <- make_panel_medz_w_cadd(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")
  )
```
```{r}
ggsave('Supplementary_revision.CADDSV_VS_WATERSHEDSV.pdf',p,height = 6,width = 6)
```




# noncoding sv performance vs cadd-sv
```{r}
# further process caddsv scores
cadd_score_pre_merge <- gtex_caddsv_score_gene_level %>%  mutate(ID = paste(SubjectID, GeneName, sep = ":")) %>% dplyr::select(ID,CADDSV_gene_score)
cadd_score_pre_merge_fibroblast <- gtex_caddsv_score_gene_level_fibroblast %>%  mutate(ID = paste(SubjectID, GeneName, sep = ":")) %>% dplyr::select(ID,CADDSV_gene_score)
cadd_score_pre_merge_muscle <- gtex_caddsv_score_gene_level_muscle %>%  mutate(ID = paste(SubjectID, GeneName, sep = ":")) %>% dplyr::select(ID,CADDSV_gene_score)
cadd_score_pre_merge_blood <- gtex_caddsv_score_gene_level_blood %>%  mutate(ID = paste(SubjectID, GeneName, sep = ":")) %>% dplyr::select(ID,CADDSV_gene_score)

# step 1: merge caddsv into watershed object. 
medz_obj <- readRDS(watershed_ABC_100k_obj_Median_file)
medz_features <- medz_obj$feat_test %>% as.data.frame() %>% rownames_to_column(var = "ID")
medz_obj$n2_posteriors$caddsv_n2_score <- medz_features %>% dplyr::select(ID) %>% left_join(cadd_score_pre_merge, by = "ID") %>% dplyr::pull(CADDSV_gene_score) %>% matrix(ncol = 1)

fibroblast_obj <- readRDS(watershed_ABC_100k_obj_Cells_Cultured_Fibroblasts_file)
fibroblast_features <- fibroblast_obj$feat_test %>% as.data.frame() %>% rownames_to_column(var = "ID")
fibroblast_obj$n2_posteriors$caddsv_n2_score <- fibroblast_features %>% dplyr::select(ID) %>% left_join(cadd_score_pre_merge_fibroblast, by = "ID") %>% dplyr::pull(CADDSV_gene_score) %>% matrix(ncol = 1)

muscle_obj <- readRDS(watershed_ABC_100k_obj_Muscle_Skeletal_file)
muscle_features <- muscle_obj$feat_test %>% as.data.frame() %>% rownames_to_column(var = "ID")
muscle_obj$n2_posteriors$caddsv_n2_score <- muscle_features %>% dplyr::select(ID) %>% left_join(cadd_score_pre_merge_muscle, by = "ID") %>% dplyr::pull(CADDSV_gene_score) %>% matrix(ncol = 1)

blood_obj <- readRDS(watershed_ABC_100k_obj_Whole_Blood_file)
blood_features <- blood_obj$feat_test %>% as.data.frame() %>% rownames_to_column(var = "ID")
blood_obj$n2_posteriors$caddsv_n2_score <- blood_features %>% dplyr::select(ID) %>% left_join(cadd_score_pre_merge_blood, by = "ID") %>% dplyr::pull(CADDSV_gene_score) %>% matrix(ncol = 1)

watershed_ABC_100k_obj_Cells_Cultured_Fibroblasts_noncoding <- re_calculate_precision_recall_caddsv_noncoding(fibroblast_obj,data_Cells_Cultured_Fibroblasts_file,1,0.2,0.1,compute_ci = FALSE)
watershed_ABC_100k_obj_Muscle_Skeletal_noncoding <- re_calculate_precision_recall_caddsv_noncoding(muscle_obj,data_Muscle_Skeletal_file,1,0.2,0.1,compute_ci = FALSE)
watershed_ABC_100k_obj_Whole_Blood_noncoding <- re_calculate_precision_recall_caddsv_noncoding(blood_obj,data_Whole_Blood_file,1,0.2,0.1,compute_ci = FALSE)
watershed_ABC_100k_obj_MedZ_noncoding <- re_calculate_precision_recall_caddsv_noncoding(medz_obj,data_Median_file,1,0.2,0.1,compute_ci = FALSE)

# watershed
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') %>% arrange(desc(row_number()))
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') %>% arrange(desc(row_number()))
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') %>% arrange(desc(row_number()))
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') %>% arrange(desc(row_number()))
# gam
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 = 'WGS-only',outlier_type = 'Blood') %>% arrange(desc(row_number()))
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 = 'WGS-only',outlier_type = 'Fibroblast') %>% arrange(desc(row_number()))
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 = 'WGS-only',outlier_type = 'Muscle') %>% arrange(desc(row_number()))
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 = 'WGS-only',outlier_type = 'MedZ') %>% arrange(desc(row_number()))
# caddsv
blood_caddsv_noncoding <- data_frame(recall = watershed_ABC_100k_obj_Whole_Blood_noncoding[[1]]$evaROC$CADDSV_recall,
           precision = watershed_ABC_100k_obj_Whole_Blood_noncoding[[1]]$evaROC$CADDSV_precision,
           prediction_type = 'CADDSV',outlier_type = 'Blood') %>% arrange(desc(row_number()))
fibroblast_caddsv_noncoding <- data_frame(recall = watershed_ABC_100k_obj_Cells_Cultured_Fibroblasts_noncoding[[1]]$evaROC$CADDSV_recall,
           precision = watershed_ABC_100k_obj_Cells_Cultured_Fibroblasts_noncoding[[1]]$evaROC$CADDSV_precision,
           prediction_type = 'CADDSV',outlier_type = 'Fibroblast') %>% arrange(desc(row_number()))
muscle_caddsv_noncoding <- data_frame(recall = watershed_ABC_100k_obj_Muscle_Skeletal_noncoding[[1]]$evaROC$CADDSV_recall,
           precision = watershed_ABC_100k_obj_Muscle_Skeletal_noncoding[[1]]$evaROC$CADDSV_precision,
           prediction_type = 'CADDSV',outlier_type = 'Muscle') %>% arrange(desc(row_number()))
medz_caddsv_noncoding <- data_frame(recall = watershed_ABC_100k_obj_MedZ_noncoding[[1]]$evaROC$CADDSV_recall,
           precision = watershed_ABC_100k_obj_MedZ_noncoding[[1]]$evaROC$CADDSV_precision,
           prediction_type = 'CADDSV',outlier_type = 'MedZ') %>% arrange(desc(row_number()))

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,
                                 fibroblast_caddsv_noncoding, blood_caddsv_noncoding, muscle_caddsv_noncoding, medz_caddsv_noncoding)

p_noncoding <- make_panel_4d_caddsv(total_pr_data_noncoding) + labs('Noncoding performance in different outlier signals') + theme(plot.margin = margin(2,2,2,2,"cm"))
```
```{r}
p_noncoding
```



# doing PR-curve by distance
```{r}
# further process caddsv scores
cadd_score_pre_merge <- gtex_caddsv_score_gene_level %>%  mutate(ID = paste(SubjectID, GeneName, sep = ":")) %>% dplyr::select(ID,CADDSV_gene_score)
cadd_score_pre_merge_fibroblast <- gtex_caddsv_score_gene_level_fibroblast %>%  mutate(ID = paste(SubjectID, GeneName, sep = ":")) %>% dplyr::select(ID,CADDSV_gene_score)
cadd_score_pre_merge_muscle <- gtex_caddsv_score_gene_level_muscle %>%  mutate(ID = paste(SubjectID, GeneName, sep = ":")) %>% dplyr::select(ID,CADDSV_gene_score)
cadd_score_pre_merge_blood <- gtex_caddsv_score_gene_level_blood %>%  mutate(ID = paste(SubjectID, GeneName, sep = ":")) %>% dplyr::select(ID,CADDSV_gene_score)

# step 1: merge caddsv into watershed object. 
medz_obj <- readRDS(watershed_ABC_100k_obj_Median_file)
medz_features <- medz_obj$feat_test %>% as.data.frame() %>% rownames_to_column(var = "ID")
medz_obj$n2_posteriors$caddsv_n2_score <- medz_features %>% dplyr::select(ID) %>% left_join(cadd_score_pre_merge, by = "ID") %>% dplyr::pull(CADDSV_gene_score) %>% matrix(ncol = 1)

fibroblast_obj <- readRDS(watershed_ABC_100k_obj_Cells_Cultured_Fibroblasts_file)
fibroblast_features <- fibroblast_obj$feat_test %>% as.data.frame() %>% rownames_to_column(var = "ID")
fibroblast_obj$n2_posteriors$caddsv_n2_score <- fibroblast_features %>% dplyr::select(ID) %>% left_join(cadd_score_pre_merge_fibroblast, by = "ID") %>% dplyr::pull(CADDSV_gene_score) %>% matrix(ncol = 1)

muscle_obj <- readRDS(watershed_ABC_100k_obj_Muscle_Skeletal_file)
muscle_features <- muscle_obj$feat_test %>% as.data.frame() %>% rownames_to_column(var = "ID")
muscle_obj$n2_posteriors$caddsv_n2_score <- muscle_features %>% dplyr::select(ID) %>% left_join(cadd_score_pre_merge_muscle, by = "ID") %>% dplyr::pull(CADDSV_gene_score) %>% matrix(ncol = 1)

blood_obj <- readRDS(watershed_ABC_100k_obj_Whole_Blood_file)
blood_features <- blood_obj$feat_test %>% as.data.frame() %>% rownames_to_column(var = "ID")
blood_obj$n2_posteriors$caddsv_n2_score <- blood_features %>% dplyr::select(ID) %>% left_join(cadd_score_pre_merge_blood, by = "ID") %>% dplyr::pull(CADDSV_gene_score) %>% matrix(ncol = 1)
```

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

    # Make Precision-recall curve for ASE
    # This plot contains 3 curves (one for Watershed, RIVER, and GAM)
	distance_category <- 'intragenic coding or noncoding'
	outlier_name <- distance_category
  	plotter_intragenic <- ggplot(data=df[as.character(df$distance_category)==distance_category,], 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","#D81B60")) +
                scale_linetype_manual(values=c("solid", "solid","solid")) +
                gtex_v8_figure_theme() + 
                draw_label(outlier_name,x=.5,y=.95,size=12) + 
  	            theme(plot.title = element_text(size=20))

    
    # Make Precision-recall curve for Expression
    # This plot contains 3 curves (one for Watershed, RIVER, and GAM)
	distance_category <- '1b-50kb noncoding'
	outlier_name <- distance_category
  	plotter_1_10k <- ggplot(data=df[as.character(df$distance_category)==distance_category,], 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","#D81B60")) +
                scale_linetype_manual(values=c("solid", "solid","solid")) +
                gtex_v8_figure_theme() +
                draw_label(outlier_name,x=.5,y=.95,size=12) + 
  	            theme(plot.title = element_text(size=20))
  	# Make Precision-recall curve for Expression
    # This plot contains 3 curves (one for Watershed, RIVER, and GAM)
	distance_category <- '50kb-100kb noncoding'
	outlier_name <- distance_category
  	plotter_10k_100k <- ggplot(data=df[as.character(df$distance_category)==distance_category,], 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","#D81B60")) +
                scale_linetype_manual(values=c("solid", "solid","solid")) +
                gtex_v8_figure_theme() + 
                draw_label(outlier_name,x=.5,y=.95,size=12) + 
  	            theme(plot.title = element_text(size=20))


    # Get legend from ASE Plot (Note: this is the same legend as found in the the Expression and Splicing plots)
    legend <- get_legend(plotter_intragenic)

    # Combine PR-curves from Expression, ASE, and splicing into 1 plot via cowplot
    combined_plots <- plot_grid(plotter_intragenic + theme(legend.position="none"), 
                                plotter_1_10k + theme(legend.position="none"), 
                                plotter_10k_100k+ theme(legend.position="bottom"),
                                rel_heights = c(1,1,1), ncol = 1)
    #combined_w_legend <- plot_grid(combined_plots, legend, ncol = 1, rel_heights = c(0.95, 0.05))


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

	return(combined)
}
```


```{r}
watershed_ABC_100k_obj_Cells_Cultured_Fibroblasts_dist_bin_1 <- re_calculate_precision_recall_caddsv_by_distance(fibroblast_obj,data_Cells_Cultured_Fibroblasts_file,1,0.2,0.1,compute_ci = FALSE, distance_bin = 'intragenic coding or noncoding')
watershed_ABC_100k_obj_Cells_Cultured_Fibroblasts_dist_bin_2 <- re_calculate_precision_recall_caddsv_by_distance(fibroblast_obj,data_Cells_Cultured_Fibroblasts_file,1,0.2,0.1,compute_ci = FALSE, distance_bin = '1b-50kb noncoding')
watershed_ABC_100k_obj_Cells_Cultured_Fibroblasts_dist_bin_3 <- re_calculate_precision_recall_caddsv_by_distance(fibroblast_obj,data_Cells_Cultured_Fibroblasts_file,1,0.2,0.1,compute_ci = FALSE, distance_bin = '50kb-100kb noncoding')

```
```{r}
# watershed
fibroblast_watershed_1 <- data_frame(recall = watershed_ABC_100k_obj_Cells_Cultured_Fibroblasts_dist_bin_1[[1]]$evaROC$watershed_recall,
           precision = watershed_ABC_100k_obj_Cells_Cultured_Fibroblasts_dist_bin_1[[1]]$evaROC$watershed_precision,
           prediction_type = 'Watershed', distance_category = 'intragenic coding or noncoding',outlier_type = 'Fibroblast') %>% arrange(desc(row_number()))
fibroblast_watershed_2 <- data_frame(recall = watershed_ABC_100k_obj_Cells_Cultured_Fibroblasts_dist_bin_2[[1]]$evaROC$watershed_recall,
           precision = watershed_ABC_100k_obj_Cells_Cultured_Fibroblasts_dist_bin_2[[1]]$evaROC$watershed_precision,
           prediction_type = 'Watershed', distance_category = '1b-50kb noncoding',outlier_type = 'Fibroblast') %>% arrange(desc(row_number()))
fibroblast_watershed_3 <- data_frame(recall = watershed_ABC_100k_obj_Cells_Cultured_Fibroblasts_dist_bin_3[[1]]$evaROC$watershed_recall,
           precision = watershed_ABC_100k_obj_Cells_Cultured_Fibroblasts_dist_bin_3[[1]]$evaROC$watershed_precision,
           prediction_type = 'Watershed', distance_category = '50kb-100kb noncoding',outlier_type = 'Fibroblast') %>% arrange(desc(row_number()))

# gam
fibroblast_gam_1 <- data_frame(recall = watershed_ABC_100k_obj_Cells_Cultured_Fibroblasts_dist_bin_1[[1]]$evaROC$GAM_recall,
           precision = watershed_ABC_100k_obj_Cells_Cultured_Fibroblasts_dist_bin_1[[1]]$evaROC$GAM_precision,
           prediction_type = 'WGS-only', distance_category = 'intragenic coding or noncoding',outlier_type = 'Fibroblast') %>% arrange(desc(row_number()))
fibroblast_gam_2 <- data_frame(recall = watershed_ABC_100k_obj_Cells_Cultured_Fibroblasts_dist_bin_2[[1]]$evaROC$GAM_recall,
           precision = watershed_ABC_100k_obj_Cells_Cultured_Fibroblasts_dist_bin_2[[1]]$evaROC$GAM_precision,
           prediction_type = 'WGS-only', distance_category = '1b-50kb noncoding',outlier_type = 'Fibroblast') %>% arrange(desc(row_number()))
fibroblast_gam_3 <- data_frame(recall = watershed_ABC_100k_obj_Cells_Cultured_Fibroblasts_dist_bin_3[[1]]$evaROC$GAM_recall,
           precision = watershed_ABC_100k_obj_Cells_Cultured_Fibroblasts_dist_bin_3[[1]]$evaROC$GAM_precision,
           prediction_type = 'WGS-only', distance_category = '50kb-100kb noncoding',outlier_type = 'Fibroblast') %>% arrange(desc(row_number()))

# caddsv
fibroblast_caddsv_1 <- data_frame(recall = watershed_ABC_100k_obj_Cells_Cultured_Fibroblasts_dist_bin_1[[1]]$evaROC$CADDSV_recall,
           precision = watershed_ABC_100k_obj_Cells_Cultured_Fibroblasts_dist_bin_1[[1]]$evaROC$CADDSV_precision,
           prediction_type = 'CADDSV', distance_category = 'intragenic coding or noncoding',outlier_type = 'Fibroblast') %>% arrange(desc(row_number()))
fibroblast_caddsv_2 <- data_frame(recall = watershed_ABC_100k_obj_Cells_Cultured_Fibroblasts_dist_bin_2[[1]]$evaROC$CADDSV_recall,
           precision = watershed_ABC_100k_obj_Cells_Cultured_Fibroblasts_dist_bin_2[[1]]$evaROC$CADDSV_precision,
           prediction_type = 'CADDSV', distance_category = '1b-50kb noncoding',outlier_type = 'Fibroblast') %>% arrange(desc(row_number()))
fibroblast_caddsv_3 <- data_frame(recall = watershed_ABC_100k_obj_Cells_Cultured_Fibroblasts_dist_bin_3[[1]]$evaROC$CADDSV_recall,
           precision = watershed_ABC_100k_obj_Cells_Cultured_Fibroblasts_dist_bin_3[[1]]$evaROC$CADDSV_precision,
           prediction_type = 'CADDSV', distance_category = '50kb-100kb noncoding',outlier_type = 'Fibroblast') %>% arrange(desc(row_number()))


fibroblast_pr_by_distance <- rbind(fibroblast_watershed_1,fibroblast_watershed_2,fibroblast_watershed_3,
                                   fibroblast_gam_1,fibroblast_gam_2,fibroblast_gam_3,
                                   fibroblast_caddsv_1,fibroblast_caddsv_2,fibroblast_caddsv_3)

p_distance_fibro <- make_panel_4d_caddsv_distance(fibroblast_pr_by_distance) + labs('fibroblast performance comparison by distance bin') + theme(plot.margin = margin(2,2,2,2,"cm"))
```





```{r}
ggsave('Supplementary_revision.fibroblast_watershed_by_distance.png',p_distance_fibro,height = 5,width = 15)
```

```{r}
watershed_ABC_100k_obj_Cells_Cultured_Fibroblasts_noncoding <- re_calculate_precision_recall_caddsv_by_distance(fibroblast_obj,data_Cells_Cultured_Fibroblasts_file,1,0.2,0.1,compute_ci = FALSE)
watershed_ABC_100k_obj_Muscle_Skeletal_noncoding <- re_calculate_precision_recall_caddsv_by_distance(muscle_obj,data_Muscle_Skeletal_file,1,0.2,0.1,compute_ci = FALSE,)
watershed_ABC_100k_obj_Whole_Blood_noncoding <- re_calculate_precision_recall_caddsv_by_distance(blood_obj,data_Whole_Blood_file,1,0.2,0.1,compute_ci = FALSE)
watershed_ABC_100k_obj_MedZ_noncoding <- re_calculate_precision_recall_caddsv_by_distance(medz_obj,data_Median_file,1,0.2,0.1,compute_ci = FALSE)

# watershed
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') %>% arrange(desc(row_number()))
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') %>% arrange(desc(row_number()))
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') %>% arrange(desc(row_number()))
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') %>% arrange(desc(row_number()))
# gam
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 = 'WGS-only',outlier_type = 'Blood') %>% arrange(desc(row_number()))
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 = 'WGS-only',outlier_type = 'Fibroblast') %>% arrange(desc(row_number()))
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 = 'WGS-only',outlier_type = 'Muscle') %>% arrange(desc(row_number()))
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 = 'WGS-only',outlier_type = 'MedZ') %>% arrange(desc(row_number()))
# caddsv
blood_caddsv_noncoding <- data_frame(recall = watershed_ABC_100k_obj_Whole_Blood_noncoding[[1]]$evaROC$CADDSV_recall,
           precision = watershed_ABC_100k_obj_Whole_Blood_noncoding[[1]]$evaROC$CADDSV_precision,
           prediction_type = 'CADDSV',outlier_type = 'Blood') %>% arrange(desc(row_number()))
fibroblast_caddsv_noncoding <- data_frame(recall = watershed_ABC_100k_obj_Cells_Cultured_Fibroblasts_noncoding[[1]]$evaROC$CADDSV_recall,
           precision = watershed_ABC_100k_obj_Cells_Cultured_Fibroblasts_noncoding[[1]]$evaROC$CADDSV_precision,
           prediction_type = 'CADDSV',outlier_type = 'Fibroblast') %>% arrange(desc(row_number()))
muscle_caddsv_noncoding <- data_frame(recall = watershed_ABC_100k_obj_Muscle_Skeletal_noncoding[[1]]$evaROC$CADDSV_recall,
           precision = watershed_ABC_100k_obj_Muscle_Skeletal_noncoding[[1]]$evaROC$CADDSV_precision,
           prediction_type = 'CADDSV',outlier_type = 'Muscle') %>% arrange(desc(row_number()))
medz_caddsv_noncoding <- data_frame(recall = watershed_ABC_100k_obj_MedZ_noncoding[[1]]$evaROC$CADDSV_recall,
           precision = watershed_ABC_100k_obj_MedZ_noncoding[[1]]$evaROC$CADDSV_precision,
           prediction_type = 'CADDSV',outlier_type = 'MedZ') %>% arrange(desc(row_number()))

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,
                                 fibroblast_caddsv_noncoding, blood_caddsv_noncoding, muscle_caddsv_noncoding, medz_caddsv_noncoding)

p_noncoding <- make_panel_4d_caddsv(total_pr_data_noncoding) + labs('Noncoding performance in different outlier signals') + theme(plot.margin = margin(2,2,2,2,"cm"))
```

```{r}
watershed_ABC_100k_obj_Whole_Blood_dist_bin_1 <- re_calculate_precision_recall_caddsv_by_distance(blood_obj,data_Whole_Blood_file,1,0.2,0.05,compute_ci = FALSE, distance_bin = 'intragenic coding or noncoding')
watershed_ABC_100k_obj_Whole_Blood_dist_bin_2 <- re_calculate_precision_recall_caddsv_by_distance(blood_obj,data_Whole_Blood_file,1,0.2,0.05,compute_ci = FALSE, distance_bin = '1b-50kb noncoding')
watershed_ABC_100k_obj_Whole_Blood_dist_bin_3 <- re_calculate_precision_recall_caddsv_by_distance(blood_obj,data_Whole_Blood_file,1,0.2,0.05,compute_ci = FALSE, distance_bin = '50kb-100kb noncoding')

```
```{r}
# watershed
blood_watershed_1 <- data_frame(recall = watershed_ABC_100k_obj_Whole_Blood_dist_bin_1[[1]]$evaROC$watershed_recall,
           precision = watershed_ABC_100k_obj_Whole_Blood_dist_bin_1[[1]]$evaROC$watershed_precision,
           prediction_type = 'Watershed', distance_category = 'intragenic coding or noncoding',outlier_type = 'Blood') %>% arrange(desc(row_number()))
blood_watershed_2 <- data_frame(recall = watershed_ABC_100k_obj_Whole_Blood_dist_bin_2[[1]]$evaROC$watershed_recall,
           precision = watershed_ABC_100k_obj_Whole_Blood_dist_bin_2[[1]]$evaROC$watershed_precision,
           prediction_type = 'Watershed', distance_category = '1b-50kb noncoding',outlier_type = 'Blood') %>% arrange(desc(row_number()))
blood_watershed_3 <- data_frame(recall = watershed_ABC_100k_obj_Whole_Blood_dist_bin_3[[1]]$evaROC$watershed_recall,
           precision = watershed_ABC_100k_obj_Whole_Blood_dist_bin_3[[1]]$evaROC$watershed_precision,
           prediction_type = 'Watershed', distance_category = '50kb-100kb noncoding',outlier_type = 'Blood') %>% arrange(desc(row_number()))

# gam
blood_gam_1 <- data_frame(recall = watershed_ABC_100k_obj_Whole_Blood_dist_bin_1[[1]]$evaROC$GAM_recall,
           precision = watershed_ABC_100k_obj_Whole_Blood_dist_bin_1[[1]]$evaROC$GAM_precision,
           prediction_type = 'WGS-only', distance_category = 'intragenic coding or noncoding',outlier_type = 'Blood') %>% arrange(desc(row_number()))
blood_gam_2 <- data_frame(recall = watershed_ABC_100k_obj_Whole_Blood_dist_bin_2[[1]]$evaROC$GAM_recall,
           precision = watershed_ABC_100k_obj_Whole_Blood_dist_bin_2[[1]]$evaROC$GAM_precision,
           prediction_type = 'WGS-only', distance_category = '1b-50kb noncoding',outlier_type = 'Blood') %>% arrange(desc(row_number()))
blood_gam_3 <- data_frame(recall = watershed_ABC_100k_obj_Whole_Blood_dist_bin_3[[1]]$evaROC$GAM_recall,
           precision = watershed_ABC_100k_obj_Whole_Blood_dist_bin_3[[1]]$evaROC$GAM_precision,
           prediction_type = 'WGS-only', distance_category = '50kb-100kb noncoding',outlier_type = 'Blood') %>% arrange(desc(row_number()))

# caddsv
blood_caddsv_1 <- data_frame(recall = watershed_ABC_100k_obj_Whole_Blood_dist_bin_1[[1]]$evaROC$CADDSV_recall,
           precision = watershed_ABC_100k_obj_Whole_Blood_dist_bin_1[[1]]$evaROC$CADDSV_precision,
           prediction_type = 'CADDSV', distance_category = 'intragenic coding or noncoding',outlier_type = 'Blood') %>% arrange(desc(row_number()))
blood_caddsv_2 <- data_frame(recall = watershed_ABC_100k_obj_Whole_Blood_dist_bin_2[[1]]$evaROC$CADDSV_recall,
           precision = watershed_ABC_100k_obj_Whole_Blood_dist_bin_2[[1]]$evaROC$CADDSV_precision,
           prediction_type = 'CADDSV', distance_category = '1b-50kb noncoding',outlier_type = 'Blood') %>% arrange(desc(row_number()))
blood_caddsv_3 <- data_frame(recall = watershed_ABC_100k_obj_Whole_Blood_dist_bin_3[[1]]$evaROC$CADDSV_recall,
           precision = watershed_ABC_100k_obj_Whole_Blood_dist_bin_3[[1]]$evaROC$CADDSV_precision,
           prediction_type = 'CADDSV', distance_category = '50kb-100kb noncoding',outlier_type = 'Blood') %>% arrange(desc(row_number()))


blood_pr_by_distance <- rbind(blood_watershed_1,blood_watershed_2,blood_watershed_3,
                                   blood_gam_1,blood_gam_2,blood_gam_3,
                                   blood_caddsv_1,blood_caddsv_2,blood_caddsv_3)

p_distance_blood <- make_panel_4d_caddsv_distance(blood_pr_by_distance) + labs('blood performance comparison by distance bin') + theme(plot.margin = margin(2,2,2,2,"cm"))
```


```{r}
ggsave('Supplementary_revision.blood_watershed_by_distance.png',p_distance_blood,height = 5,width = 15)
```

```{r}
watershed_ABC_100k_obj_Cells_Cultured_Fibroblasts_dist_bin_1 <- re_calculate_precision_recall_caddsv_by_distance(blood_obj,data_Whole_Blood_file,1,0.2,0.05,compute_ci = FALSE, distance_bin = 'intragenic coding or noncoding')
watershed_ABC_100k_obj_Whole_Blood_dist_bin_2 <- re_calculate_precision_recall_caddsv_by_distance(blood_obj,data_Whole_Blood_file,1,0.2,0.05,compute_ci = FALSE, distance_bin = '1b-50kb noncoding')
watershed_ABC_100k_obj_Whole_Blood_dist_bin_3 <- re_calculate_precision_recall_caddsv_by_distance(blood_obj,data_Whole_Blood_file,1,0.2,0.05,compute_ci = FALSE, distance_bin = '50kb-100kb noncoding')

```
```{r}
# watershed
blood_watershed_1 <- data_frame(recall = watershed_ABC_100k_obj_Whole_Blood_dist_bin_1[[1]]$evaROC$watershed_recall,
           precision = watershed_ABC_100k_obj_Whole_Blood_dist_bin_1[[1]]$evaROC$watershed_precision,
           prediction_type = 'Watershed', distance_category = 'intragenic coding or noncoding',outlier_type = 'Blood') %>% arrange(desc(row_number()))
blood_watershed_2 <- data_frame(recall = watershed_ABC_100k_obj_Whole_Blood_dist_bin_2[[1]]$evaROC$watershed_recall,
           precision = watershed_ABC_100k_obj_Whole_Blood_dist_bin_2[[1]]$evaROC$watershed_precision,
           prediction_type = 'Watershed', distance_category = '1b-50kb noncoding',outlier_type = 'Blood') %>% arrange(desc(row_number()))
blood_watershed_3 <- data_frame(recall = watershed_ABC_100k_obj_Whole_Blood_dist_bin_3[[1]]$evaROC$watershed_recall,
           precision = watershed_ABC_100k_obj_Whole_Blood_dist_bin_3[[1]]$evaROC$watershed_precision,
           prediction_type = 'Watershed', distance_category = '50kb-100kb noncoding',outlier_type = 'Blood') %>% arrange(desc(row_number()))

# gam
blood_gam_1 <- data_frame(recall = watershed_ABC_100k_obj_Whole_Blood_dist_bin_1[[1]]$evaROC$GAM_recall,
           precision = watershed_ABC_100k_obj_Whole_Blood_dist_bin_1[[1]]$evaROC$GAM_precision,
           prediction_type = 'WGS-only', distance_category = 'intragenic coding or noncoding',outlier_type = 'Blood') %>% arrange(desc(row_number()))
blood_gam_2 <- data_frame(recall = watershed_ABC_100k_obj_Whole_Blood_dist_bin_2[[1]]$evaROC$GAM_recall,
           precision = watershed_ABC_100k_obj_Whole_Blood_dist_bin_2[[1]]$evaROC$GAM_precision,
           prediction_type = 'WGS-only', distance_category = '1b-50kb noncoding',outlier_type = 'Blood') %>% arrange(desc(row_number()))
blood_gam_3 <- data_frame(recall = watershed_ABC_100k_obj_Whole_Blood_dist_bin_3[[1]]$evaROC$GAM_recall,
           precision = watershed_ABC_100k_obj_Whole_Blood_dist_bin_3[[1]]$evaROC$GAM_precision,
           prediction_type = 'WGS-only', distance_category = '50kb-100kb noncoding',outlier_type = 'Blood') %>% arrange(desc(row_number()))

# caddsv
blood_caddsv_1 <- data_frame(recall = watershed_ABC_100k_obj_Whole_Blood_dist_bin_1[[1]]$evaROC$CADDSV_recall,
           precision = watershed_ABC_100k_obj_Whole_Blood_dist_bin_1[[1]]$evaROC$CADDSV_precision,
           prediction_type = 'CADDSV', distance_category = 'intragenic coding or noncoding',outlier_type = 'Blood') %>% arrange(desc(row_number()))
blood_caddsv_2 <- data_frame(recall = watershed_ABC_100k_obj_Whole_Blood_dist_bin_2[[1]]$evaROC$CADDSV_recall,
           precision = watershed_ABC_100k_obj_Whole_Blood_dist_bin_2[[1]]$evaROC$CADDSV_precision,
           prediction_type = 'CADDSV', distance_category = '1b-50kb noncoding',outlier_type = 'Blood') %>% arrange(desc(row_number()))
blood_caddsv_3 <- data_frame(recall = watershed_ABC_100k_obj_Whole_Blood_dist_bin_3[[1]]$evaROC$CADDSV_recall,
           precision = watershed_ABC_100k_obj_Whole_Blood_dist_bin_3[[1]]$evaROC$CADDSV_precision,
           prediction_type = 'CADDSV', distance_category = '50kb-100kb noncoding',outlier_type = 'Blood') %>% arrange(desc(row_number()))


blood_pr_by_distance <- rbind(blood_watershed_1,blood_watershed_2,blood_watershed_3,
                                   blood_gam_1,blood_gam_2,blood_gam_3,
                                   blood_caddsv_1,blood_caddsv_2,blood_caddsv_3)

p_distance_blood <- make_panel_4d_caddsv_distance(blood_pr_by_distance) + labs('blood performance comparison by distance bin') + theme(plot.margin = margin(2,2,2,2,"cm"))
```


```{r}
watershed_ABC_100k_obj_MedZ_dist_bin_1 <- re_calculate_precision_recall_caddsv_by_distance(medz_obj,data_Median_file,1,0.176,0.05,compute_ci = FALSE, distance_bin = 'intragenic coding or noncoding')
watershed_ABC_100k_obj_MedZ_dist_bin_2 <- re_calculate_precision_recall_caddsv_by_distance(medz_obj,data_Median_file,1,0.176,0.05,compute_ci = FALSE, distance_bin = '1b-50kb noncoding')
watershed_ABC_100k_obj_MedZ_dist_bin_3 <- re_calculate_precision_recall_caddsv_by_distance(medz_obj,data_Median_file,1,0.176,0.05,compute_ci = FALSE, distance_bin = '50kb-100kb noncoding')

```


```{r}
# watershed
med_watershed_1 <- data_frame(recall = watershed_ABC_100k_obj_MedZ_dist_bin_1[[1]]$evaROC$watershed_recall,
           precision = watershed_ABC_100k_obj_MedZ_dist_bin_1[[1]]$evaROC$watershed_precision,
           prediction_type = 'Watershed', distance_category = 'intragenic coding or noncoding',outlier_type = 'MedZ') %>% arrange(desc(row_number()))
med_watershed_2 <- data_frame(recall = watershed_ABC_100k_obj_MedZ_dist_bin_2[[1]]$evaROC$watershed_recall,
           precision = watershed_ABC_100k_obj_MedZ_dist_bin_2[[1]]$evaROC$watershed_precision,
           prediction_type = 'Watershed', distance_category = '1b-50kb noncoding',outlier_type = 'MedZ') %>% arrange(desc(row_number()))
med_watershed_3 <- data_frame(recall = watershed_ABC_100k_obj_MedZ_dist_bin_3[[1]]$evaROC$watershed_recall,
           precision = watershed_ABC_100k_obj_MedZ_dist_bin_3[[1]]$evaROC$watershed_precision,
           prediction_type = 'Watershed', distance_category = '50kb-100kb noncoding',outlier_type = 'MedZ') %>% arrange(desc(row_number()))

# gam
med_gam_1 <- data_frame(recall = watershed_ABC_100k_obj_MedZ_dist_bin_1[[1]]$evaROC$GAM_recall,
           precision = watershed_ABC_100k_obj_MedZ_dist_bin_1[[1]]$evaROC$GAM_precision,
           prediction_type = 'WGS-only', distance_category = 'intragenic coding or noncoding',outlier_type = 'MedZ') %>% arrange(desc(row_number()))
med_gam_2 <- data_frame(recall = watershed_ABC_100k_obj_MedZ_dist_bin_2[[1]]$evaROC$GAM_recall,
           precision = watershed_ABC_100k_obj_MedZ_dist_bin_2[[1]]$evaROC$GAM_precision,
           prediction_type = 'WGS-only', distance_category = '1b-50kb noncoding',outlier_type = 'MedZ') %>% arrange(desc(row_number()))
med_gam_3 <- data_frame(recall = watershed_ABC_100k_obj_MedZ_dist_bin_3[[1]]$evaROC$GAM_recall,
           precision = watershed_ABC_100k_obj_MedZ_dist_bin_3[[1]]$evaROC$GAM_precision,
           prediction_type = 'WGS-only', distance_category = '50kb-100kb noncoding',outlier_type = 'MedZ') %>% arrange(desc(row_number()))

# caddsv
med_caddsv_1 <- data_frame(recall = watershed_ABC_100k_obj_MedZ_dist_bin_1[[1]]$evaROC$CADDSV_recall,
           precision = watershed_ABC_100k_obj_MedZ_dist_bin_1[[1]]$evaROC$CADDSV_precision,
           prediction_type = 'CADDSV', distance_category = 'intragenic coding or noncoding',outlier_type = 'MedZ') %>% arrange(desc(row_number()))
med_caddsv_2 <- data_frame(recall = watershed_ABC_100k_obj_MedZ_dist_bin_2[[1]]$evaROC$CADDSV_recall,
           precision = watershed_ABC_100k_obj_MedZ_dist_bin_2[[1]]$evaROC$CADDSV_precision,
           prediction_type = 'CADDSV', distance_category = '1b-50kb noncoding',outlier_type = 'MedZ') %>% arrange(desc(row_number()))
med_caddsv_3 <- data_frame(recall = watershed_ABC_100k_obj_MedZ_dist_bin_3[[1]]$evaROC$CADDSV_recall,
           precision = watershed_ABC_100k_obj_MedZ_dist_bin_3[[1]]$evaROC$CADDSV_precision,
           prediction_type = 'CADDSV', distance_category = '50kb-100kb noncoding',outlier_type = 'MedZ') %>% arrange(desc(row_number()))


medz_pr_by_distance <- rbind(med_watershed_1,med_watershed_2,med_watershed_3,
                                   med_gam_1,med_gam_2,med_gam_3,
                                   med_caddsv_1,med_caddsv_2,med_caddsv_3)

p_distance_medz <- make_panel_4d_caddsv_distance(medz_pr_by_distance) + 
  labs('Medz performance comparison by distance bin') + 
  theme(plot.margin = margin(2,2,2,2,"cm"))
theme_set(theme_cowplot(font_size=5))
```



```{r}
p_distance_medz
```


```{r}
ggsave('Supplementary_revision.medz_watershed_by_distance.large_text.png',p_distance_medz,height = 12,width = 6)
```
