## ----setup, include=FALSE-----------------------------------------------------
library(PscanR)

knitr::opts_chunk$set(
    collapse = TRUE,
    comment = "#>",
    fig.align = "center",
    warning = FALSE,
    message = FALSE
)

run <- readRDS(system.file(
    "extdata", "vignettes", "custom_background", "ucbg_summary.rds",
    package = "PscanR"
))
verification <- run$verification
fullbg <- run$fullbg

timing_of <- function(manifest, step) {
    manifest$elapsed_seconds[manifest$step == step]
}
timings_1core <- read.delim(
    system.file(
        "extdata", "vignettes", "custom_background", "timings_1core.tsv",
        package = "PscanR"
    ),
    stringsAsFactors = FALSE
)
timings_8core <- read.delim(
    system.file(
        "extdata", "vignettes", "custom_background", "timings_8core.tsv",
        package = "PscanR"
    ),
    stringsAsFactors = FALSE
)
build_speedup <- round(
    timing_of(timings_1core, "background_short") /
        timing_of(timings_8core, "background_short"), 1
)
full_bg_all_mb <- round(
    fullbg$size_bytes * verification$motif_count /
        length(fullbg$motifs) / 1024^2
)
is_cerevisiae <- vapply(
    TFBSTools::getMatrixSet(
        JASPAR2020::JASPAR2020,
        list(collection = "CORE", tax_group = "fungi")
    ),
    function(m) any(grepl("cerevisiae", TFBSTools::tags(m)$species)),
    logical(1)
)
n_cerevisiae <- sum(is_cerevisiae)
pct_cerevisiae <- round(100 * n_cerevisiae / length(is_cerevisiae))

## ----installation, eval=FALSE-------------------------------------------------
# BiocManager::install("PscanR")

## ----available----------------------------------------------------------------
head(ps_available_bg(), 4)

## ----small-background-roundtrip-----------------------------------------------
example_sequences <- readRDS(system.file(
    "extdata", "vignettes", "custom_background", "rp_promoters.rds",
    package = "PscanR"
))
n_example <- min(40L, length(example_sequences))
example_sequences <- example_sequences[seq_len(n_example)]
example_motifs <- JASPAR2020::JASPAR2020
example_motifs <- TFBSTools::getMatrixSet(
    example_motifs, list(collection = "CORE", tax_group = "fungi")
)[1:2]
short_example <- ps_build_bg(example_sequences, example_motifs)
short_file <- tempfile(fileext = ".txt")
ps_write_bg_to_file(short_example, short_file)
short_restored <- ps_retrieve_bg_from_file(short_file, example_motifs)
stopifnot(isTRUE(all.equal(ps_get_bg_table(short_example),
            ps_get_bg_table(short_restored))))
full_example <- ps_build_bg(example_sequences, example_motifs, fullBG = TRUE)
full_file <- tempfile(fileext = ".rds")
saveRDS(full_example, full_file, compress = "xz")
full_restored <- readRDS(full_file)
stopifnot(identical(full_example, full_restored))
query <- names(example_sequences)[1:4]
queried <- pscan_full_bg(query, full_restored)
ps_results_table(queried)
unlink(c(short_file, full_file))

## ----promoters, eval=FALSE----------------------------------------------------
# 
# # UCSC's sgdGene table, packaged and versioned.
# library(TxDb.Scerevisiae.UCSC.sacCer3.sgdGene)
# txdb <- TxDb.Scerevisiae.UCSC.sacCer3.sgdGene
# 
# # The background covers the 16 nuclear chromosomes.
# tx <- GenomicFeatures::transcripts(txdb, columns = c("tx_id", "tx_name"))
# nuclear <- paste0("chr", as.roman(seq_len(16)))
# tx <- tx[as.character(GenomicRanges::seqnames(tx)) %in% nuclear]
# names(tx) <- GenomicRanges::mcols(tx)$tx_name
# 
# # promoters() is strand-aware and anchored on the transcription start site.
# ranges <- GenomicRanges::promoters(tx, upstream = 950, downstream = 50)
# 
# genome <- BSgenome.Scerevisiae.UCSC.sacCer3::BSgenome.Scerevisiae.UCSC.sacCer3
# 
# # Exclude out-of-bounds ranges before sequence extraction; getSeq() otherwise
# # raises an error.
# lengths <- GenomeInfoDb::seqlengths(genome)[
#     as.character(GenomicRanges::seqnames(ranges))
# ]
# inside <- GenomicRanges::start(ranges) >= 1L &
#     GenomicRanges::end(ranges) <= lengths
# ranges <- ranges[inside]
# 
# prom_seq <- Biostrings::getSeq(genome, ranges)
# names(prom_seq) <- names(ranges)

## ----matrices-----------------------------------------------------------------
motifs <- TFBSTools::getMatrixSet(
    JASPAR2020::JASPAR2020,
    list(collection = "CORE", tax_group = "fungi")
)
length(motifs)

## ----build, eval=FALSE--------------------------------------------------------
# 
# background <- ps_build_bg(
#     prom_seq, motifs,
#     BPPARAM = BiocParallel::SerialParam()
# )

## ----timing-table, echo=FALSE-------------------------------------------------
build_timings <- data.frame(
    Step = sprintf("ps_build_bg(), %d matrices", verification$motif_count),
    `One core` = timing_of(timings_1core, "background_short"),
    `Eight cores` = timing_of(timings_8core, "background_short"),
    check.names = FALSE
)
build_timings$Speedup <- round(
    build_timings$`One core` / build_timings$`Eight cores`, 1
)
knitr::kable(
    build_timings, digits = 1,
    caption = paste(
        "Background construction time in seconds,",
        "using one or eight cores."
    )
)

## ----write, eval=FALSE--------------------------------------------------------
# ps_write_bg_to_file(background, "J2020_sacCer3_950u_50d_UCSC.rebuild.txt")

## ----show-file----------------------------------------------------------------
bg_file <- system.file(
    "extdata", "vignettes", "custom_background",
    "J2020_sacCer3_950u_50d_UCSC.rebuild.txt",
    package = "PscanR"
)
cat(readLines(bg_file, n = 4), sep = "\n")
file.size(bg_file)

## ----read-back----------------------------------------------------------------
background <- ps_retrieve_bg_from_file(bg_file, motifs)
background
head(ps_get_bg_table(background), 3)

## ----verify-table-------------------------------------------------------------
knitr::kable(
    data.frame(
        Check = c(
            "Matrices compared",
            "BG_SIZE equal",
            "BG_MEAN agreeing to 10 decimal places",
            "BG_STDEV agreeing to 10 decimal places"
        ),
        Result = c(
            verification$motif_count,
            verification$size_equal,
            sprintf("%d / %d", verification$mean_agree_10dp,
                verification$motif_count),
            sprintf("%d / %d", verification$stdev_agree_10dp,
                verification$motif_count)
        )
    ),
    align = c("l", "r")
)

## ----build-full, eval=FALSE---------------------------------------------------
# 
# # The bundled fullbg_motifs.tsv records the fifteen selected matrices.
# subset_ids <- read.delim(system.file(
#     "extdata", "vignettes", "custom_background", "fullbg_motifs.tsv",
#     package = "PscanR"
# ))$MATRIX
# 
# full_bg <- ps_build_bg(
#     prom_seq, motifs[subset_ids],
#     BPPARAM = BiocParallel::SerialParam(), fullBG = TRUE
# )

## ----save-full, eval=FALSE----------------------------------------------------
# saveRDS(full_bg, "sacCer3_950u_50d_fullBG_15motifs.rds", compress = "xz")

## ----load-full----------------------------------------------------------------
full_bg <- readRDS(system.file(
    "extdata", "vignettes", "custom_background",
    "sacCer3_950u_50d_fullBG_15motifs.rds",
    package = "PscanR"
))
TFBSTools::name(full_bg)

## ----foreground---------------------------------------------------------------
rp_promoters <- readRDS(system.file(
    "extdata", "vignettes", "custom_background", "rp_promoters.rds",
    package = "PscanR"
))
rp_ids <- names(rp_promoters)
length(rp_ids)

## ----run-fullbg, message=TRUE-------------------------------------------------
by_legend <- pscan_full_bg(rp_ids, full_bg)

## ----fullbg-results-----------------------------------------------------------
enrichment <- ps_results_table(by_legend)
head(enrichment[, c("NAME", "ZSCORE", "P.VALUE", "FDR")], 5)

## ----workflow, eval=FALSE-----------------------------------------------------
# 
# # 1. Extract promoters using the intended foreground window.
# 
# genome <- BSgenome.Scerevisiae.UCSC.sacCer3::BSgenome.Scerevisiae.UCSC.sacCer3
# 
# library(TxDb.Scerevisiae.UCSC.sacCer3.sgdGene)
# txdb <- TxDb.Scerevisiae.UCSC.sacCer3.sgdGene
# 
# tx <- GenomicFeatures::transcripts(txdb, columns = c("tx_id", "tx_name"))
# nuclear <- paste0("chr", as.roman(seq_len(16)))
# tx <- tx[as.character(GenomicRanges::seqnames(tx)) %in% nuclear]
# names(tx) <- GenomicRanges::mcols(tx)$tx_name
# 
# ranges <- GenomicRanges::promoters(tx, upstream = 950, downstream = 50)
# 
# lengths <- GenomeInfoDb::seqlengths(genome)[
#     as.character(GenomicRanges::seqnames(ranges))
# ]
# 
# inside <- GenomicRanges::start(ranges) >= 1L &
#     GenomicRanges::end(ranges) <= lengths
# 
# ranges <- ranges[inside]
# 
# prom_seq <- Biostrings::getSeq(genome, ranges)
# 
# # 2. Select motif matrices for the organism's taxonomic group.
# motifs <- TFBSTools::getMatrixSet(
#     JASPAR2020::JASPAR2020, list(collection = "CORE", tax_group = "fungi")
# )
# 
# # 3. Construct the standard background.
# background <- ps_build_bg(
#     prom_seq, motifs, BPPARAM = BiocParallel::SerialParam()
# )
# 
# # 4. Save the background statistics.
# ps_write_bg_to_file(background, "my_background.txt")
# background <- ps_retrieve_bg_from_file("my_background.txt", motifs)
# 
# # Optionally, retain all hits in a full background for identifier-based queries.
# full_bg <- ps_build_bg(prom_seq, motifs, fullBG = TRUE)
# saveRDS(full_bg, "my_full_background.rds", compress = "xz")
# full_bg <- readRDS("my_full_background.rds")
# 
# # Query identifiers must match the background's promoter naming scheme.
# # This example uses SGD systematic identifiers for cytoplasmic ribosomal
# # protein genes. An external list can be read from a file with readLines().
# gene_names <- unlist(as.list(org.Sc.sgd.db::org.Sc.sgdGENENAME))
# my_gene_ids <- names(gene_names)[grepl("^RP[LS][0-9]", gene_names)]
# 
# ps_results_table(pscan_full_bg(my_gene_ids, full_bg))

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