## ----quick-start--------------------------------------------------------------
# exons <- prepare_exons(
#   txdb = <A TxDb OBJECT>,
#   dtu_table = <DTU_TABLE>,
#   coef_col = "estimate"
# )
# 
# exons <- preprocess(exons, coef_col = "estimate")
# 
# # find skipped exons:
# skipped <- exons |> find_se()

## ----dtu-table----------------------------------------------------------------
library(readr)
# load DTU results
dir <- system.file("extdata", package="splicelogic")
dtu_table <- readr::read_delim(file.path(dir, "dtu_table.tsv"))

## ----load-exons---------------------------------------------------------------
library(plyranges)
exons_file <- "exons_M31.bed.gz"
exons_mcols_file <- "exons_mcols.tsv.gz"
exons <- plyranges::read_bed(file.path(dir, exons_file))
mcols(exons) <- DataFrame(
  readr::read_delim(file.path(dir, exons_mcols_file))
)

## ----add-seqinfo--------------------------------------------------------------
si <- Seqinfo::Seqinfo(genome="mm39")
seqlevels(exons) <- seqlevels(si)
seqinfo(exons) <- si

## ----merge-dtu----------------------------------------------------------------
txp_idx <- match(exons$tx_id, dtu_table$tx_id)
cols_to_add <- dtu_table[txp_idx, ] |> dplyr::select(-tx_id)
merged_DF <- S4Vectors::cbind(mcols(exons), cols_to_add)
mcols(exons) <- merged_DF

## ----join-mcols---------------------------------------------------------------
# exons <- exons |>
#   plyranges::join_mcols_left(dtu_table, by = "tx_id")

## ----filter-sig---------------------------------------------------------------
sig_exons <- exons |> filter(padj < .1)

## ----preprocess---------------------------------------------------------------
library(splicelogic)
sig_exons <- sig_exons |>
  preprocess(coef_col = "estimate")

## ----calc-se------------------------------------------------------------------
skipped <- sig_exons |> find_se()
skipped

## ----calc-re------------------------------------------------------------------
included <- sig_exons |> find_ie()
included

## ----calc-mx------------------------------------------------------------------
mx <- sig_exons |> find_mxe()
mx

## ----calc-ri------------------------------------------------------------------
ri <- sig_exons |> find_ri()
ri

## ----calc-a5-a3-ss------------------------------------------------------------
a5ss <- sig_exons |> find_a5ss()
a5ss
a3ss <- sig_exons |> find_a3ss()
a3ss

## ----calc-all-----------------------------------------------------------------
all_events <- sig_exons |> find_all_events()
all_events

## ----events-barplot-----------------------------------------------------------
barplot(
  table(all_events$event), horiz=TRUE, las=1,
  xlab="exons participating in an event"
)

## ----ahub---------------------------------------------------------------------
suppressPackageStartupMessages({
  library(AnnotationHub)
  library(AnnotationDbi)
  library(GenomicFeatures)
})
ah <- AnnotationHub()
txdb <- ah[["AH75191"]] # GENCODE v32 (human)

## ----get-txps-----------------------------------------------------------------
suppressPackageStartupMessages({
  library(tibble)
})
txps <- txdb |>
  AnnotationDbi::select(
    keys(txdb, "TXID"), c("TXNAME","GENEID"), "TXID"
  ) |>
  tibble::as_tibble() |>
  dplyr::select(tx_num = TXID, tx_id = TXNAME, gene_id = GENEID) |>
  dplyr::filter(!is.na(gene_id))
txps

## ----sim-dtu------------------------------------------------------------------
# simulate DTU results
sim_dtu_table <- txps |>
  dplyr::mutate(
    padj = runif(dplyr::n()),
    effect_est = rnorm(dplyr::n())
  )
sim_dtu_table

## ----prepare-exons------------------------------------------------------------
human_exons <- prepare_exons(
  txdb, sim_dtu_table, coef_col = "effect_est", verbose = TRUE
)
human_exons <- human_exons |>
  filter(padj < .01) |>
  preprocess(coef_col = "effect_est")

## ----extract-exons------------------------------------------------------------
# extract exons as a GRangesList
exons_list <- GenomicFeatures::exonsBy(
  txdb,
  by="tx"
  )
# Our DTU table aligns with txps, which aligns with the names
# of the GRangesList. prepare_exons() handles alignment checks.
names(exons_list) <- sim_dtu_table$tx_id

## ----flatten-exons------------------------------------------------------------
flat_exons <- unlist(exons_list)
# swap tx_id with exon_name as the names of the GRanges
flat_exons$tx_id <- names(flat_exons) # store transcript ids
names(flat_exons) <- flat_exons$exon_name

## ----add-dtu------------------------------------------------------------------
txp_idx <- match(flat_exons$tx_id, sim_dtu_table$tx_id)
cols_to_add <- sim_dtu_table[txp_idx,] |>
  dplyr::select(-c(tx_id, tx_num))
merged_DF <- cbind(mcols(flat_exons), cols_to_add)
mcols(flat_exons) <- merged_DF

## ----session-info-------------------------------------------------------------
sessionInfo()