Add SpliceAI API, comprehensive tests, fix ggplot2 facet_grid

- New SPLICEAI_query_api(): query Broad SpliceAI Lookup REST API
  for on-the-fly splice predictions without local precomputed files
- Fix deprecated facet_grid(facets=) -> rows/cols/vars() across all
  IMPACT and DEEPLEARNING plot functions
- Fix IMPACT_get_annotations CHR type mismatch (char vs int)
- Fix IMPACT_iterate_get_annotations to handle errors gracefully
- Fix DEEPLEARNING_plot combn() crash with < 2 SNP groups
- Fix DEEPLEARNING_plot outlier removal edge case
- 57 tests passing, covering all 25 exported functions with real
  data integration tests (Zenodo, GitHub, SpliceAI API)

Co-Authored-By: Claude Opus 4.6 (1M context) <noreply@anthropic.com>

Author: bschilder

DESCRIPTION

@@ -59,7 +59,8 @@ Suggests:
     shades,
     reshape2,
     stringr,
-    R.utils
+    R.utils,
+    jsonlite
 Remotes:
     github::RajLabMSSM/echodata,
     github::RajLabMSSM/echotabix,

NAMESPACE

@@ -19,6 +19,7 @@ export(IMPACT_postprocess_annotations)
 export(IMPACT_query)
 export(IMPACT_snp_group_boxplot)
 export(SPLICEAI_plot)
+export(SPLICEAI_query_api)
 export(SPLICEAI_query_tsv)
 export(SPLICEAI_query_tsv_iterate)
 export(SPLICEAI_query_vcf)

R/DEEPLEARNING_plot.R

@@ -72,7 +72,10 @@ DEEPLEARNING_plot <- function(
         ## Suppress the side-effect plot device from boxplot()
         outliers <- graphics::boxplot(annot_melt$value,
                                       plot = FALSE)$out
-        annot_melt <- annot_melt[-which(annot_melt$value %in% outliers), ]
+        if (length(outliers) > 0) {
+            idx <- which(annot_melt$value %in% outliers)
+            if (length(idx) > 0) annot_melt <- annot_melt[-idx, ]
+        }
     }
     colorDict <- echodata::snp_group_colorDict()
     dat_plot <- subset(
@@ -84,14 +87,19 @@ DEEPLEARNING_plot <- function(
                                levels = names(colorDict),
                                ordered = TRUE)
         )
-    snp_groups_present <- unique(dat_plot$SNP_group)
-    comparisons <- utils::combn(
-        x = as.character(snp_groups_present),
-        m = 2,
-        FUN = comparisons_filter,
-        simplify = FALSE
-    )
-    comparisons <- Filter(Negate(is.null), comparisons)
+    snp_groups_present <- unique(as.character(
+        dat_plot$SNP_group[!is.na(dat_plot$SNP_group)]))
+    if (length(snp_groups_present) >= 2) {
+        comparisons <- utils::combn(
+            x = snp_groups_present,
+            m = 2,
+            FUN = comparisons_filter,
+            simplify = FALSE
+        )
+        comparisons <- Filter(Negate(is.null), comparisons)
+    } else {
+        comparisons <- list()
+    }
 
     method <- "wilcox.test"
     gp <- ggplot2::ggplot(
@@ -102,7 +110,7 @@ DEEPLEARNING_plot <- function(
         ggplot2::geom_violin(alpha = 0.6) +
         ggplot2::geom_boxplot(alpha = 0.6, color = "black") +
         ggplot2::facet_grid(
-            facets = stats::as.formula(facet_formula),
+            rows = stats::as.formula(facet_formula),
             scales = "free"
         ) +
         ggplot2::labs(

R/IMPACT_get_annotations.R

@@ -39,7 +39,7 @@ IMPACT_get_annotations <- function(
     Annot <- IMPACT_score <- POS <- NULL;
 
     if (!is.null(dat)) {
-        dat$CHR <- gsub("chr", "", dat$CHR)
+        dat$CHR <- as.integer(gsub("chr", "", dat$CHR))
         chrom <- dat$CHR[1]
     }
     URL <- file.path(baseURL,
@@ -126,7 +126,7 @@ IMPACT_iterate_get_annotations <- function(
         function(chrom, .merged_DT = merged_DT) {
             messager("+ IMPACT:: Gathering annotations for chrom =",
                      chrom, v = verbose)
-            try({
+            tryCatch({
                 dat <- subset(.merged_DT, CHR == chrom)
                 annot_melt <- IMPACT_get_annotations(
                     baseURL = baseURL,
@@ -153,8 +153,13 @@ IMPACT_iterate_get_annotations <- function(
                          "annotations found at IMPACT_score >=",
                          IMPACT_score_thresh, v = verbose)
                 return(annot_melt)
+            }, error = function(e) {
+                messager("+ IMPACT:: Error for chrom", chrom, ":",
+                         conditionMessage(e), v = verbose)
+                return(NULL)
             })
         })
-    ANNOT_MELT <- data.table::rbindlist(ANNOT_MELT)
+    ANNOT_MELT <- Filter(Negate(is.null), ANNOT_MELT)
+    ANNOT_MELT <- data.table::rbindlist(ANNOT_MELT, fill = TRUE)
     return(ANNOT_MELT)
 }

R/IMPACT_plot.R

@@ -182,7 +182,7 @@ IMPACT_plot_enrichment <- function(ENRICH,
                              height = 0) +
         ggplot2::geom_hline(yintercept = 1, linetype = "dashed",
                             alpha = 0.8) +
-        ggplot2::facet_grid(facets = SNP.group ~ .,
+        ggplot2::facet_grid(rows = ggplot2::vars(SNP.group),
                             switch = "y", space = "free_x",
                             scales = "free_x") +
         ggplot2::theme_bw() +
@@ -196,7 +196,7 @@ IMPACT_plot_enrichment <- function(ENRICH,
         ggplot2::geom_col(position = "dodge") +
         ggplot2::geom_hline(yintercept = 1, linetype = "dashed",
                             alpha = 0.8) +
-        ggplot2::facet_grid(facets = . ~ Tissue,
+        ggplot2::facet_grid(cols = ggplot2::vars(Tissue),
                             switch = "y", space = "free_x",
                             scales = "free_x") +
         ggplot2::theme_bw() +
@@ -298,7 +298,7 @@ IMPACT_plot_impact_score <- function(annot_melt,
         subset(annot_melt, IMPACT_score > 0.5),
         ggplot2::aes(x = Mb, y = IMPACT_score, color = TF)) +
         ggplot2::geom_point(show.legend = FALSE) +
-        ggplot2::facet_grid(facets = Tissue ~ ., switch = "y") +
+        ggplot2::facet_grid(rows = ggplot2::vars(Tissue), switch = "y") +
         ggplot2::theme_bw() +
         ggplot2::labs(y = "IMPACT score per tissue") +
         ggplot2::theme(

R/SPLICEAI_query_api.R

@@ -0,0 +1,95 @@
+#' Query SpliceAI Lookup API
+#'
+#' Query the public
+#' \href{https://spliceailookup.broadinstitute.org/}{SpliceAI Lookup}
+#' web API hosted by the Broad Institute. This service runs the SpliceAI
+#' model on-the-fly for individual variants without requiring local
+#' installation or precomputed files.
+#'
+#' @param variants Character vector of variants in
+#'   \code{"chrom-pos-ref-alt"} format (e.g. \code{"2-179415988-C-CA"}).
+#' @param genome Genome build: \code{"37"} (GRCh37/hg19) or
+#'   \code{"38"} (GRCh38/hg38).
+#' @param distance Maximum distance between the variant and gained/lost
+#'   splice site (default 50, max 10000).
+#' @param mask Mask scores for annotated splice sites
+#'   (0 = raw, 1 = masked). Default 0.
+#' @param verbose Print messages.
+#'
+#' @returns A \code{data.table} with one row per variant, containing
+#'   columns for \code{variant}, \code{gene}, \code{DS_AG}, \code{DS_AL},
+#'   \code{DS_DG}, \code{DS_DL} (delta scores), and corresponding
+#'   \code{DP_*} (delta positions). Returns \code{NULL} for variants
+#'   with no scores.
+#'
+#' @details
+#' \strong{Rate limit:} The public API supports only a handful of
+#' queries per user per minute. For batch processing, use
+#' \code{\link{SPLICEAI_query_tsv}} with local precomputed files or
+#' set up a local instance via
+#' \url{https://github.com/broadinstitute/SpliceAI-lookup}.
+#'
+#' @export
+#' @family SPLICEAI
+#' @source
+#' \url{https://spliceailookup.broadinstitute.org/}
+#' \doi{10.1016/j.cell.2018.12.015}
+#' @importFrom data.table data.table rbindlist
+#' @examples
+#' \dontrun{
+#' result <- SPLICEAI_query_api(
+#'     variants = "2-179415988-C-CA",
+#'     genome = "37"
+#' )
+#' }
+SPLICEAI_query_api <- function(variants,
+                               genome = c("37", "38"),
+                               distance = 50,
+                               mask = 0,
+                               verbose = TRUE) {
+
+    genome <- match.arg(genome)
+    base_url <- paste0("https://spliceai-", genome,
+                       "-xwkwwwxdwq-uc.a.run.app/spliceai/")
+
+    results <- lapply(variants, function(v) {
+        messager("SPLICEAI:: Querying API for:", v, v = verbose)
+        url <- paste0(base_url,
+                      "?hg=", genome,
+                      "&distance=", distance,
+                      "&mask=", mask,
+                      "&variant=", v)
+        tryCatch({
+            resp <- readLines(url(url), warn = FALSE)
+            json <- jsonlite::fromJSON(resp, flatten = TRUE)
+            if (!is.null(json$error)) {
+                messager("SPLICEAI:: No scores for", v, ":", json$error,
+                         v = verbose)
+                return(NULL)
+            }
+            scores <- json$scores
+            if (is.null(scores) || length(scores) == 0) return(NULL)
+            dt <- data.table::data.table(
+                variant = v,
+                gene = scores$g_name,
+                transcript = scores$t_id,
+                DS_AG = as.numeric(scores$DS_AG),
+                DS_AL = as.numeric(scores$DS_AL),
+                DS_DG = as.numeric(scores$DS_DG),
+                DS_DL = as.numeric(scores$DS_DL),
+                DP_AG = as.integer(scores$DP_AG),
+                DP_AL = as.integer(scores$DP_AL),
+                DP_DG = as.integer(scores$DP_DG),
+                DP_DL = as.integer(scores$DP_DL)
+            )
+            return(dt)
+        }, error = function(e) {
+            messager("SPLICEAI:: API error for", v, ":",
+                     conditionMessage(e), v = verbose)
+            return(NULL)
+        })
+    })
+    results <- Filter(Negate(is.null), results)
+    if (length(results) == 0) return(NULL)
+    data.table::rbindlist(results, fill = TRUE)
+}