---
title: "Building a background for PscanR"
author:
    - Federico Zambelli
date: "`r Sys.Date()`"
output:
    BiocStyle::html_document:
        toc: true
        toc_depth: 2
vignette: >
    %\VignetteIndexEntry{5. Building a background for PscanR}
    %\VignetteEncoding{UTF-8}
    %\VignetteEngine{knitr::rmarkdown}
editor_options:
    markdown:
        wrap: 80
---

```{r setup, include=FALSE}
library(PscanR)

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

fixture <- function(...) {
    system.file(
        "extdata", "vignettes", "custom_background", ...,
        package = "PscanR", mustWork = TRUE
    )
}

run <- readRDS(fixture("ucbg_summary.rds"))
verification <- run$verification
fullbg <- run$fullbg

timing_of <- function(manifest, step) {
    manifest$elapsed_seconds[manifest$step == step]
}
timings_1core <- read.delim(
    fixture("timings_1core.tsv"), stringsAsFactors = FALSE
)
timings_8core <- read.delim(
    fixture("timings_8core.tsv"), 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

PscanR is being prepared for Bioconductor submission. After acceptance, install
it with the following command (not executed while building this vignette):

```{r installation, eval=FALSE}
BiocManager::install("PscanR")
```

For pre-submission installation from source, see the package README.

# Introduction

PscanR provides access to precomputed backgrounds for five model organisms:
human, mouse, and fruit fly (RefSeq annotation), *Arabidopsis* (TAIR
annotation),
and yeast (SGD annotation).

`ps_available_bg()` lists the available resources:

```{r available}
head(ps_available_bg(), 4)
```

This vignette describes custom background construction. A custom background is
required when the catalogue does not include the desired:

- organism or reference annotation;
- promoter window;
- motif collection.

The example reconstructs an existing catalogue entry for *Saccharomyces
cerevisiae*, using promoters from −950 to +50 bp relative to the transcription
start site (TSS). This permits comparison with the published background.

The bundled timing and verification results were generated using the developers'
computational environment. The code below documents the construction procedure; 
expensive or externally dependent steps are not evaluated during vignette 
compilation.

The yeast example permits a relatively short demonstration. In the recorded
run, background construction took
`r round(timing_of(timings_1core, "background_short"))` seconds on one core, for
`r format(verification$promoter_count, big.mark = ",")` promoters and
`r verification$motif_count` matrices. Runtime for other organisms depends on
the number and length of promoters, the motif collection, and the hardware.

# A small executable construction and round trip

The following code builds both background representations from a subset of
real yeast promoters and two JASPAR motifs. This restricted background is an
API demonstration, not a replacement for the complete promoter universe used
for enrichment. Both computations run during vignette builds.

```{r small-background-roundtrip}
example_sequences <- readRDS(fixture("rp_promoters.rds"))
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))
```

# Part A: Building a standard background

## The promoters

A background is constructed from reference sequences and motif matrices. Here,
the reference comprises annotated nuclear promoters, defined using the same
window as the intended foreground: **−950 to +50 bp relative to the TSS**.

The reference is specified by three components:

1. **Annotation:** UCSC's `sgdGene` table, provided by the
    `TxDb.Scerevisiae.UCSC.sacCer3.sgdGene` package, whose version records
    the annotation release used;
2. **Genome:** `BSgenome.Scerevisiae.UCSC.sacCer3`, whose chromosome names
    must agree with the annotation's;
3. **Window:** 950 bp upstream and 50 bp downstream of the TSS.

Changes to the assembly or annotation can alter promoter sequences and
background statistics. Record the reference versions and use a consistent
promoter definition for foreground and background.

```{r 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)
```

This produces `r format(verification$transcript_count, big.mark = ",")`
promoter sequences, one per retained transcript.

Retain the original sequence names and duplicate sequences when passing the
extracted promoters to `ps_build_bg()`. The function deduplicates sequences,
retains those of the maximum supplied length, and removes those with more
than 50% ambiguous `N` bases. Supply the intended full-length promoter window;
the function does not infer it from genomic coordinates.
Preserving names before deduplication is particularly important for full
backgrounds, where the identifier mapping must retain equivalent transcripts.
Out-of-bounds genomic ranges must still be excluded before sequence extraction,
as shown above.

The background size is the number of distinct promoter sequences remaining
after these operations: `r format(verification$promoter_count, big.mark = ",")`
in this example.

## The matrices

Supply the motif collection as a `PFMatrixList`, as defined by TFBSTools. The
catalogued backgrounds use JASPAR CORE matrices restricted to the organism's
taxonomic group.

```{r matrices}
motifs <- TFBSTools::getMatrixSet(
    JASPAR2020::JASPAR2020,
    list(collection = "CORE", tax_group = "fungi")
)
length(motifs)
```

## Building a standard background

`ps_build_bg()` retains the best-scoring window for each promoter–motif pair,
considering both strands. It then computes the background mean and standard
deviation of these scores and records the number of retained promoters.

```{r build, eval=FALSE}

background <- ps_build_bg(
    prom_seq, motifs,
    BPPARAM = BiocParallel::SerialParam()
)
```

```{r 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."
    )
)
```

Using eight cores (`BPPARAM` set to `BiocParallel::MulticoreParam(workers = 8)`)
produced a `r build_speedup`-fold speedup for the full motif collection in the
recorded run. Parallelism is distributed over motifs. Actual speedup also
depends on promoter count and length, worker overhead, and available resources.

## Saving and reloading a standard background

```{r write, eval=FALSE}
ps_write_bg_to_file(background, "J2020_sacCer3_950u_50d_UCSC.rebuild.txt")
```

The following output shows the header, the first three motif records, and the
file size in bytes:

```{r show-file}
bg_file <- fixture("J2020_sacCer3_950u_50d_UCSC.rebuild.txt")
cat(readLines(bg_file, n = 4), sep = "\n")
file.size(bg_file)
```

The text format contains a header followed by one row per matrix, with four
fields: matrix identifier, promoter count, mean score, and score standard
deviation. Storing summary statistics keeps standard background files small.

`ps_retrieve_bg_from_file()` associates these statistics with the supplied
matrices by identifier and returns a `PSMatrixList` for use with `pscan()`.

```{r read-back}
background <- ps_retrieve_bg_from_file(bg_file, motifs)
background
head(ps_get_bg_table(background), 3)
```

`ps_get_bg_table()` retrieves the three statistical columns from the reloaded
object. Text serialisation preserves their values to the precision of the
file format; it does not guarantee bit-for-bit equality with the original
in-memory floating-point values.

## Verification against the catalogue

The reconstructed and catalogued backgrounds use the same annotation, genome,
promoter window, and matrices. The following table reports their agreement in
the recorded verification run.


```{r 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")
)
```

All matrices agree in background size, and both score statistics agree to ten
decimal places in this comparison.

Floating-point arithmetic and text serialisation can introduce small numerical
differences. Agreement to ten decimal places is the criterion used for this
verification, not a guarantee of exact byte-equality across platforms. Those
tiny 
differences across platforms are usually negligible and do not impact PscanR 
results in any meaningful way. 

# Part B: Building and querying a full background

A standard background stores three summary statistics per matrix. Foreground
analysis therefore requires promoter sequences extracted using the matching
reference and window.

PscanR also supports full backgrounds, built with `fullBG = TRUE`. These retain
the per-promoter scan results, allowing subsequent foreground analyses to
retrieve existing hits by identifier. This avoids repeated DNA scanning and is
useful when analysing many gene sets against the same background.

A full background retains the best-hit score, position, strand, and matched
oligonucleotide for every retained promoter and matrix. Its identifier legend
maps original sequence names to the representatives retained after
deduplication.

`pscan_full_bg()` accepts transcript identifiers matching this legend, retrieves
the stored hits, and calculates foreground enrichment without rescanning DNA.

```{r build-full, eval=FALSE}

# The bundled fullbg_motifs.tsv records the fifteen selected matrices.
subset_ids <- read.delim(fixture("fullbg_motifs.tsv"))$MATRIX

full_bg <- ps_build_bg(
    prom_seq, motifs[subset_ids],
    BPPARAM = BiocParallel::SerialParam(), fullBG = TRUE
)
```


**Pass sequences before deduplication.** The legend is constructed from the
input names. Deduplicating beforehand prevents the function from retaining
identifiers associated with the removed duplicate sequences.

**Storage requirements.** Full backgrounds can use sizeable storage up to
hundreds 
of megabytes depending on the number of promoters and matrices. 
The bundled yeast demonstration contains `r length(fullbg$motifs)`
matrices and occupies approximately `r round(fullbg$size_bytes / 1024)` KiB.
Linear
extrapolation by matrix count gives approximately `r full_bg_all_mb` MiB for
all `r verification$motif_count` yeast matrices. Complete full-background 
collections are not maintaned nor bundled with the package.

## Saving and reloading a full background

Save a full background as an `.rds` file to preserve its `PSMatrixList`, hits,
and identifier legend. 

```{r save-full, eval=FALSE}
saveRDS(full_bg, "sacCer3_950u_50d_fullBG_15motifs.rds", compress = "xz")
```

The reported file sizes use `compress = "xz"`. This compression method can
reduce storage requirements at the cost of additional compression time.

Reloading restores the hits and legend used in the following foreground query.

```{r load-full}
full_bg <- readRDS(fixture("sacCer3_950u_50d_fullBG_15motifs.rds"))
TFBSTools::name(full_bg)
```

## Analysing a foreground by identifier

The foreground example uses yeast cytoplasmic ribosomal protein genes to
illustrate retrieval of a promoter set expected to contain Rap1 motif matches.

The set contains `r fullbg$rp_total` cytoplasmic ribosomal protein genes
selected from `org.Sc.sgd.db` by standard names matching `^RP[LS][0-9]` and
mapped to SGD systematic identifiers. The background contains
`r fullbg$rp_used` of these identifiers.

```{r foreground}
rp_promoters <- readRDS(fixture("rp_promoters.rds"))
rp_ids <- names(rp_promoters)
length(rp_ids)
```

```{r run-fullbg, message=TRUE}
by_legend <- pscan_full_bg(rp_ids, full_bg)
```

The message reports the identifier scheme detected in the background.

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

RAP1 is the highest-ranked motif, with a z-score of
`r round(fullbg$top_zscore, 1)`. SFP1 ranks second in this restricted motif
collection. The displayed table reports both nominal and adjusted p-values.


This workflow performs identifier lookup and retrieval of stored hits for each
query, followed by enrichment calculation. It avoids promoter extraction and
DNA rescanning. Query identifiers must match the background legend; gene
symbols cannot be substituted unless they are the identifiers used by that
background. Here, the SGD systematic identifiers name the stored promoters.


# Part C: Checklist for a new background

Record the following choices to support reproducibility:

1. **Annotation source.** A packaged `TxDb` if Bioconductor ships one for the
    organism. Otherwise a UCSC table via `txdbmaker::makeTxDbFromUCSC()`, or a
    GFF/GTF via `txdbmaker::makeTxDbFromGFF()`; with those, record *which
    release* was used and when it was retrieved.
2. **Genome sequence.**
    A `BSgenome` package if one exists; otherwise a FASTA read with
    `Biostrings::readDNAStringSet()`. Record the assembly and ensure that
    chromosome names agree with the annotation.
3. **Window.** Specify the upstream and downstream extent relative to the TSS.
    Foreground and background must share this window and compatible reference
    sequences and annotations.
4. **Motif collection, and its taxonomic subset.** Inspect what the collection
    contains and record its release. In the JASPAR CORE fungi collection used
    here, `r pct_cerevisiae`% of matrices have *S. cerevisiae* species
    annotation.
    Assess suitability for the organism under study.
5. **Identifier scheme.** Preserve meaningful identifier suffixes. The
    `scheme` argument of `pscan_full_bg()` supports `"refseq"`, `"sgd"`,
    `"tair"`, and `"generic"`; `"auto"` detects recognised conventions.
    Generic identifiers are matched without suffix removal.

# Summary

```{r 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))
```

# Data provenance

The installed `inst/scripts/fixture-provenance.md` describes the sources,
selection rules, serialization, and reproduction recipes for the bundled
example data. `inst/DATA_SOURCES.md` records attribution and reuse terms.
The large published analyses use saved results; small workflow examples are
executed during the build. Reproducing a full analysis requires the matching
annotation and promoter window, not just the same gene identifiers.

# Session information

```{r session-info}
sessionInfo()
```