---
title: "Motif enrichment analysis with PscanR: mouse retinal promoters"
author:
    - Federico Zambelli
date: "`r Sys.Date()`"
output:
    BiocStyle::html_document:
        toc: true
        toc_depth: 2
vignette: >
    %\VignetteIndexEntry{3. PscanR: mouse retinal promoters}
    %\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", "mouse_hd_retina", ...,
        package = "PscanR", mustWork = TRUE
    )
}
```

# 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 identifies transcription factor binding motifs enriched in regulatory
sequences. Promoters of co-expressed or co-regulated genes are scanned with a
collection of motif matrices, and their mean best-match scores are compared
with those of a reference promoter population. Enriched motifs suggest
candidate regulators.

PscanR provides access to precomputed backgrounds for five model organisms:
human, mouse, and fruit fly (RefSeq annotation), *Arabidopsis* (TAIR
annotation),
and yeast (SGD annotation). Custom backgrounds can also be constructed for
these or other organisms, as described in [Building a background for
PscanR](PscanR_Custom_Background.html).

This vignette demonstrates motif enrichment analysis of mouse promoters and
the functions available for examining the results.

## The data

The example uses gene expression data from a study of the R6/1 mouse model of
Huntington's disease, which profiled retina and striatum from the same animals.

> Cano-Cano F, Martín-Loro F, Gallardo-Orihuela A, González-Montelongo MC,
> Ortuño-Miquel S, Hervás-Corpión I, de la Villa P, Ramón-Marco L,
> Navarro-Calvo J, Gómez-Jaramillo L, Arroba AI, Valor LM. *Retinal dysfunction
> in Huntington's disease mouse models concurs with local gliosis and microglia
> activation.* Scientific Reports 14:4176 (2024).
> [doi:10.1038/s41598-024-54347-8](https://doi.org/10.1038/s41598-024-54347-8)

Differential expression was computed separately in the two tissues, and the
resulting genes were divided into six sets by tissue and direction of change.
The input table contains gene identifiers and their set membership.

The examples below use the bundled results to avoid network access and
repeated genome-wide sequence extraction.

```{r read-gene-sets}
membership <- read.delim(
    fixture("paper_subset_membership.tsv"),
    stringsAsFactors = FALSE
)
str(membership)
gene_sets <- split(membership$ensembl_gene_id, membership$subset)
lengths(gene_sets)
```

```{r subset-key, echo=FALSE}
subset_labels <- c(
    a = "retina-only down", b = "shared down", c = "striatum-only down",
    d = "retina-only up", e = "shared up", f = "striatum-only up"
)
knitr::kable(
    data.frame(
        set = names(subset_labels),
        definition = unname(subset_labels),
        genes = lengths(gene_sets)[names(subset_labels)],
        row.names = NULL
    ),
    caption = "The six gene sets, by tissue and direction of change."
)
```

The rest of this vignette follows set `d`, the
`r length(gene_sets$d)` genes upregulated in retina and not in striatum, and
uses set `f`, the `r length(gene_sets$f)` genes upregulated in striatum and not
in retina, as the comparison.

# The background

PscanR precomputed backgrounds are indexed by organism, genome assembly,
promoter window and motif collection. `ps_available_bg()` lists what is
available in the catalogue bundled with the package:

```{r discover-background}
catalog <- ps_available_bg("J2020_mm10_950u_50d", details = TRUE)
knitr::kable(
    catalog[c(
        "organism", "assembly", "upstream", "downstream",
        "jaspar_release", "motif_count", "promoter_count"
    )],
    row.names = FALSE,
    caption = "The matching precomputed background."
)
```

`ps_retrieve_bg()` retrieves the background statistics
and associates them with the matching JASPAR matrices in a `PSMatrixList`.
Here `source = "zenodo"` downloads the file from Zenodo directly; by default it
is retrieved through ExperimentHub.

```{r retrieve-background}
background <- ps_retrieve_bg(
    JASPAR_matrix = "Jaspar2020",
    org = "mm",
    prom_reg = c(-950, 50),
    assembly = "mm10",
    version = 2,
    source = "zenodo"
)
length(background)
```

This returns 746 matrices from the JASPAR 2020 CORE non-redundant vertebrate
collection. Each matrix includes the background promoter count and the mean
and standard deviation of its best-match scores across those promoters.

The selected window extends from 950 bp upstream to 50 bp downstream of the
transcription start site (TSS). **Foreground promoters must use the same window
as the background**, together with a compatible assembly and annotation.

# From genes to promoters

`pscan()` accepts promoter sequences as a `DNAStringSet`. The following
procedure constructs the foreground using the background's promoter definition.

When relevant transcripts are known, their corresponding promoters can be
selected directly. Gene-level differential expression results require an
additional selection rule because a gene may have multiple transcription start
sites. Selecting one promoter per gene avoids weighting genes by their number
of annotated transcripts. A deterministic rule makes this choice reproducible,
but does not establish which promoter is active under the experimental
conditions.

`ps_select_promoters()` maps a gene list to one promoter per mapped gene using
supplied promoter sequences and a gene-to-transcript annotation table.

Identifier suffixes have different meanings across annotation systems. The
suffix in the RefSeq accession `NM_000546.6` denotes a sequence version and is
removed for matching. In the TAIR identifier `AT1G01110.2`, the suffix denotes
a transcript variant and is preserved. PscanR detects and reports the scheme
used for matching.

The Ensembl gene identifiers in this experiment are mapped to RefSeq
transcripts through `org.Mm.eg.db`. RefSeq version normalisation permits
matching unversioned annotation entries to versioned promoter names. It does
not imply that transcript sequences or start sites are unchanged between
versions; the promoter reference must still match the background.

The following bundled example contains ten genes from this dataset with
multiple candidate transcripts:

```{r promoter-scheme, message=TRUE}
example <- readRDS(fixture("promoter_selection_example.rds"))

selection <- ps_select_promoters(
    genes = example$genes,
    promoter_sequences = example$promoter_sequences,
    annotation = example$annotation,
    gene_col = "ENSEMBL", transcript_col = "REFSEQ",
    mode = "representative",
    return = "mapping",
    sequence_names = "gene_transcript"
)
```

The returned mapping records the selected promoter and the selection rule:

```{r promoter-scheme-table, echo=FALSE}
knitr::kable(
    selection[, c(
        "gene", "promoter_id", "scheme", "suffix_role",
        "transcript_class", "selection_source", "n_candidates", "decided_by"
    )],
    row.names = FALSE,
    caption = "Selected promoters and the rules used for their selection."
)
```

For RefSeq the rule is to prefer an `NM_` messenger transcript over an `NR_`
non-coding one, and to break any remaining tie by accession order. Protein
accessions (`NP_`, `XP_`) are excluded before ranking.

This dataset requires promoter sequences matching the background selected
above: mm10 assembly, the UCSC `ncbiRefSeqCurated` track, 950 bp upstream to 50
bp downstream of the TSS.

The following code illustrates promoter extraction and selection. It is not
evaluated during vignette compilation because it requires additional
annotation and genome packages: `AnnotationDbi`, `Biostrings`, `GenomicRanges`,
`IRanges`, `org.Mm.eg.db`, and `BSgenome.Mmusculus.UCSC.mm10`.

`mm10_950u_50d_refseq_snapshot.rds` is shipped with this example and records the
transcript identifiers and TSS coordinates of the promoter universe the
background was computed over.


```{r select-real, eval=FALSE}

retina_up_genes <- gene_sets$d

# Transcript identifiers and start sites used to construct the background.
snapshot <- readRDS(fixture("mm10_950u_50d_refseq_snapshot.rds"))
transcripts <- snapshot$transcripts

transcript_tss <- GenomicRanges::GRanges(
    seqnames = transcripts$seqname,
    ranges = IRanges::IRanges(start = transcripts$tss, width = 1L),
    strand = transcripts$strand
)
names(transcript_tss) <- transcripts$transcript_id

# promoters() takes the region around each 5' end, so it is strand-aware and
# anchored on the transcription start site.
ranges <- GenomicRanges::promoters(
    transcript_tss,
    upstream = snapshot$upstream,
    downstream = snapshot$downstream
)

genome <- BSgenome.Mmusculus.UCSC.mm10::BSgenome.Mmusculus.UCSC.mm10
promoters_mm10 <- Biostrings::getSeq(genome, ranges)
names(promoters_mm10) <- names(ranges)

# Map Ensembl genes to RefSeq accessions. Multiple rows can occur per gene;
# ps_select_promoters() excludes protein accessions.
annotation <- AnnotationDbi::select(
    org.Mm.eg.db::org.Mm.eg.db,
    keys = retina_up_genes,
    keytype = "ENSEMBL",
    columns = c("REFSEQ", "SYMBOL")
)

retina_mapping <- ps_select_promoters(
    genes = retina_up_genes,
    promoter_sequences = promoters_mm10,
    annotation = annotation,
    gene_col = "ENSEMBL", transcript_col = "REFSEQ",
    mode = "representative",
    return = "mapping",
    sequence_names = "gene_transcript"
)

retina_up_promoters <- promoters_mm10[retina_mapping$promoter_id]
```

The mapping and scan results are bundled in a precomputed object, allowing the
remaining examples to run without repeating promoter extraction or scanning.

```{r load-precomputed}
precomputed <- readRDS(fixture("uc2_vignette_data.rds"))
retina_mapping <- precomputed$retina_mapping
```

`ps_selection_summary()` reports the numbers of input and mapped genes,
including those with multiple candidate promoters.

```{r selection-summary}
ps_selection_summary(retina_mapping)$overall
```

Nearly half the mapped genes had more than one candidate promoter. Most of
these selections required accession-based tie-breaking after transcript
ranking. In the absence of condition-specific transcript evidence, this rule
provides a reproducible representative rather than a validated active promoter.

The unmapped genes lack a matching curated RefSeq promoter in this reference:
most have no RefSeq accession in the annotation used for the analysis, three are
mitochondrial genes that carry only a protein accession and are outside the
nuclear promoter set, and one has only predicted `XM_` models, which the
curated track excludes. The mapping records these genes in its
`ps_unmapped_genes` attribute:

```{r unmapped}
head(attr(retina_mapping, "ps_unmapped_genes"), 4)
```

`selection_source` records the rule used to select each promoter:

```{r selection-by-source}
ps_selection_summary(retina_mapping)$by_source
```

# A small executable scan

This example selects the real mm10 promoters bundled above and scans two
motifs using their matched -950/+50 background statistics. It demonstrates the
complete API with a small gene set; interpret the full retinal analysis below
for biological conclusions. FDR here applies only to these two motifs.

```{r small-scan}
small_sequences <- ps_select_promoters(
    genes = example$genes, promoter_sequences = example$promoter_sequences,
    annotation = example$annotation, gene_col = "ENSEMBL",
    transcript_col = "REFSEQ", mode = "representative", return = "sequences"
)
small_background <- precomputed$results_subset_d[1:2]
small_results <- pscan(small_sequences, small_background)
ps_results_table(small_results)
```

# Running the scan

`pscan()` takes the foreground promoter sequences and the background
`PSMatrixList`, and returns the updated `PSMatrixList` carrying the
per-sequence results for every motif.

The `BPPARAM` argument configures execution through BiocParallel:

- `SerialParam()` uses the current R process and simplifies package checks and
    debugging.

- `SnowParam()` starts independent socket-based R workers; it is the portable
    parallel choice on Windows and Unix-like systems, but has greater worker
    startup and memory overhead.

- `MulticoreParam()` uses forked workers on supported Unix-like systems. It can
    reduce data-transfer overhead, but suitability depends on the execution
    environment. Forked execution is unavailable on Windows.

```{r scan-time, include=FALSE}
timings <- precomputed$timings
scan_seconds <- round(
    sum(timings$elapsed_seconds[timings$subset %in% c("d", "f")]), 0
)
```

The following scan is not evaluated during vignette compilation. In the
recorded run, scanning the two sets used here took `r scan_seconds`
seconds with `BPPARAM` set to
`BiocParallel::MulticoreParam(workers = 8)`.

```{r scan, eval=FALSE}

BPPARAM <- BiocParallel::SerialParam() # Execute in the current R process.

# Portable parallel execution with four independent R workers:
# BPPARAM <- BiocParallel::SnowParam(workers = 4, type = "SOCK")

# Forked execution on supported Unix-like systems:
# BPPARAM <- BiocParallel::MulticoreParam(workers = 4)

results <- pscan(
    retina_up_promoters, background,
    BPPARAM = BPPARAM
)
```

`ps_results_table()` converts the result object to a motif enrichment table.

```{r results-table, eval=FALSE}
enrichment <- ps_results_table(results)
```

The following objects contain the saved scan results. Subsequent examples
operate on these bundled results.

```{r load-results}
# Saved enrichment table for retina-only upregulated promoters.
enrichment <- precomputed$enrichment_d
dim(enrichment)

# A motif subset of the PSMatrixList retains the per-promoter results.
results_subset <- precomputed$results_subset_d
results_subset
```

The subset includes interferon regulatory factors (IRFs) and STATs relevant to
the original study, the highest-ranked GC-rich motifs, and a motif with a
negative enrichment z-score for comparison.

```{r subset-names}
TFBSTools::name(results_subset)
```

The enrichment table contains one row per motif, indexed by JASPAR matrix
identifier:

```{r results-table-shown}
head(enrichment[order(-enrichment$ZSCORE), ], 8)
```

GC-rich motifs, including KLF15, ZBTB14, and SP and EGR zinc fingers, occupy the
highest ranks. IRF3, the highest-ranked IRF motif, has rank
`r which(rownames(enrichment[order(-enrichment$ZSCORE), ]) == "MA1418.1")`
with a z-score of
`r sprintf("%.2f", enrichment["MA1418.1", "ZSCORE"])`.

The background parameters remain associated with the result and can be
retrieved using `ps_get_bg_table()`:

```{r bg-table}
ps_get_bg_table(results_subset)
```

`BG_SIZE` is the number of mouse promoters the background was computed over;
`BG_MEAN` and `BG_STDEV` describe how that motif scores across them.

For each motif, the foreground mean score is compared with the background
mean using an enrichment z-score. The table reports the corresponding
upper-tail p-value and its Benjamini–Hochberg adjustment across the motif
collection in the `FDR` column.

`ps_motif_barplot()` displays the highest-ranked motifs. The `n` argument
limits their number to at most 20 by default. With `group = "class"`, bar
colours indicate transcription factor structural classes extracted from the
JASPAR metadata by `ps_motif_class()`.

```{r top-motifs, fig.width=8, fig.height=5, dpi=120}

table(ps_motif_class(results_subset))

ps_motif_barplot(results_subset, group = "class")
```

Within the subset, C2H2 zinc fingers rank above the interferon factors,
consistent with the complete table. NEUROD2 illustrates a negative z-score,
indicating a foreground mean below the background mean. Its z-score of
`r sprintf("%.1f", enrichment["MA0668.1", "ZSCORE"])`
is associated with p ≈ 1 because the reported p-value is an upper-tail test
for enrichment. 

The function returns a `ggplot` object that can be customised with a title,
theme, or additional layers:

```{r barplot-composable}
ps_motif_barplot(results_subset, group = "class") +
    ggplot2::labs(title = "Retina-only upregulated promoters")
```

# Checking the result against a second gene set

The GC-rich zinc-finger motifs and interferon-related motifs both show
enrichment, with the zinc-finger motifs having higher z-scores. Their relative
statistical significance does not establish which pattern is specific to the
biological condition.

A positive enrichment z-score indicates that a motif's mean score exceeds the
background expectation. This may reflect condition-associated regulation,
shared promoter composition, or other properties of the selected gene set.

Set `f`, comprising striatum-only upregulated genes identified by the same
differential expression pipeline, provides a comparison. Examining enrichment
in both sets helps assess whether the motif patterns are shared across these
foregrounds. This comparison is descriptive; it is not a direct statistical
test of a difference between tissues.

```{r second-set-table}
enrichment_f <- precomputed$enrichment_f

data.frame(
    set = c("d", "f"),
    definition = unname(subset_labels[c("d", "f")]),
    promoters = precomputed$subset_qc$scan_eligible_sequences[
        match(c("d", "f"), precomputed$subset_qc$subset)
    ],
    nominal = c(
        sum(enrichment$P.VALUE < 0.05), sum(enrichment_f$P.VALUE < 0.05)
    ),
    fdr = c(sum(enrichment$FDR < 0.05), sum(enrichment_f$FDR < 0.05))
)
```

The following plot compares the same motifs in the striatal foreground:

```{r second-set-barplot, fig.width=8, fig.height=5, dpi=120}
ps_motif_barplot(precomputed$results_subset_f, group = "class")
```

The GC-rich zinc-finger motifs remain enriched: the KLF15 z-score changes from
`r sprintf("%.2f", enrichment["MA1513.1", "ZSCORE"])` to
`r sprintf("%.2f", enrichment_f["MA1513.1", "ZSCORE"])`.
In contrast, the IRF3 z-score changes from
`r sprintf("%.2f", enrichment["MA1418.1", "ZSCORE"])` to
`r sprintf("%.2f", enrichment_f["MA1418.1", "ZSCORE"])`,
and every IRF and STAT motif in the panel has a negative z-score in the
striatal set.

The interferon-related pattern differs between these two tissue-associated
gene sets, consistent with the study's reported retinal inflammatory response.
The zinc-finger pattern is shared. These interpretations depend on the
comparison between foregrounds rather than on the p-values of a single scan.

The shared pattern could reflect regulators active in both tissues or common
properties of the selected promoters. Distinguishing these explanations
requires additional evidence. The present comparison supports different
patterns of motif enrichment across the selected sets, without establishing
their cause.

## Comparison with the published figure

The published analysis makes the same comparison over all six sets and all
eighteen IRF and STAT matrices, as its Supplementary Figure S2C.

The six `data.frame` objects returned by `ps_results_table()` can be combined
into a long-format table for visualisation.

```{r irf-stat-data}
irf_stat <- precomputed$irf_stat
str(irf_stat[c("subset", "motif", "ZSCORE", "significance")])
```

The `motif` ordering and nominal significance thresholds follow the published
figure to facilitate comparison. These annotations differ from the adjusted
p-values reported in the `FDR` column.

```{r irf-stat-heatmap, fig.width=10, fig.height=6.7, dpi=120}
ggplot2::ggplot(
    irf_stat,
    ggplot2::aes(x = subset, y = motif, fill = ZSCORE)
) +
    ggplot2::geom_tile(color = "white", linewidth = 0.25) +
    ggplot2::geom_text(ggplot2::aes(label = significance), size = 2.5) +
    ggplot2::scale_fill_gradient2(
        low = "#2166ac", mid = "white", high = "#b2182b",
        midpoint = 0, name = "Z-score"
    ) +
    ggplot2::labs(
        x = NULL, y = NULL,
        title = "IRF and STAT motif enrichment across the six gene sets",
        subtitle = "* nominal P < 0.05; ** nominal P < 0.005"
    ) +
    ggplot2::theme_minimal(base_size = 10) +
    ggplot2::theme(
        panel.grid = ggplot2::element_blank(),
        axis.text.x = ggplot2::element_text(angle = 35, hjust = 1),
        plot.title.position = "plot"
    )
```

In the retina-only upregulated set, 17 of 18 matrices have positive z-scores
and 12 meet an FDR threshold of 0.05. Across the three downregulated sets,
the family's mean z-score is negative, most strongly in the striatum-only set.
These additional foregrounds extend the comparison across tissues and
directions of expression change; negative scores remain descriptive measures
of depletion rather than lower-tail significance tests.

```{r irf-stat-counts, echo=FALSE}
counts <- do.call(rbind, lapply(split(irf_stat, irf_stat$subset), function(x) {
    data.frame(
        set = as.character(x$subset[[1]]),
        `nominal P < 0.05` = sum(x$P.VALUE < 0.05),
        `FDR < 0.05` = sum(x$FDR < 0.05),
        `mean Z` = round(mean(x$ZSCORE), 2),
        check.names = FALSE
    )
}))
knitr::kable(
    counts, row.names = FALSE,
    caption = paste(
        "The eighteen IRF and STAT matrices, counted in each set."
    )
)
```

```{r subset-e-size, include=FALSE}
qc <- precomputed$subset_qc
n_subset_e <- qc$scan_eligible_sequences[qc$subset == "e"]
```

Set `e` contains only `r n_subset_e` promoters. Its small size limits
interpretation of the normal approximation
used for enrichment. 

# Inspecting the result object

`ps_results_table()` summarises results across promoters. The `PSMatrixList`
also retains per-promoter results, which can be retrieved using accessor
functions.

A `PSMatrixList` can be indexed by JASPAR matrix identifier; `MA1418.1` selects
the IRF3 entry:

```{r accessors}
irf3 <- results_subset[["MA1418.1"]]
TFBSTools::name(irf3)

ps_zscore(irf3)
ps_pvalue(irf3)
c(background = ps_bg_avg(irf3), foreground = ps_fg_avg(irf3))
ps_hits_size(irf3)
```

The modest difference between foreground and background means produces an
enrichment z-score of `r sprintf("%.2f", ps_zscore(irf3))`. This statistic
concerns the mean across `r ps_hits_size(irf3)` promoters; it does not measure
the effect size or binding strength in an individual promoter.

`ps_hits_table()` returns each promoter's best motif match, including its score,
position, strand, and matched oligonucleotide.

```{r hits-table}
head(ps_hits_table(irf3, pos_shift = -950L), 5)
```

The oligonucleotide is reported in the **forward orientation of the supplied
promoter sequence**, including reverse-strand matches. The motif-orientation
sequence for a reverse-strand match is its reverse complement. The IRF3 matrix
contains a tandem `GAAAC` repeat characteristic of interferon-stimulated
response elements:

```{r consensus}
pfm <- TFBSTools::Matrix(irf3)
paste(rownames(pfm)[apply(pfm, 2, which.max)], collapse = "")
```


# Motif score patterns across promoters

`ps_z_table()` returns a matrix of standardised best-hit scores, with promoters
as rows and motifs as columns. These per-promoter values differ from the
motif-level enrichment z-scores, which compare foreground and background means.

```{r z-table}
z <- ps_z_table(results_subset)
dim(z)
round(z[1:4, 1:5], 2)
```

`ps_zscore_heatmap()` displays these values to compare motif score patterns
across promoters. By default, it scales each row and clusters both axes using
correlation distances.

```{r zscore-heatmap, fig.width=9, fig.height=7, dpi=110}
ps_zscore_heatmap(
    results_subset, FDR = 1,
    show_rownames = FALSE, treeheight_row = 0
)
```

`FDR = 1` displays every supplied motif, including the negative-score
comparison motif. The plotted subset was selected from the full 746-motif scan.

The column dendrogram separates the three C2H2 zinc-finger motifs from the
remaining motifs. The row-scaled heatmap suggests that promoters with
relatively high zinc-finger scores differ from those with relatively high
interferon-related scores. This view describes heterogeneity within the
foreground that is not apparent from motif-level enrichment alone.

# Positional distribution of motif matches

`ps_hitpos_map()` plots the position of the best hit in each promoter.

```{r hitpos-map, fig.width=9, fig.height=7, dpi=110}
ps_hitpos_map(
    results_subset, FDR = 1, shift = -950,
    show_rownames = FALSE, treeheight_row = 0
)
```

`shift = -950` relabels the axis so that coordinates are relative to the
transcription start site, matching how the promoters were extracted. The
zinc-finger motifs show many hits within approximately 200 bp of the TSS,
whereas interferon-related hits are distributed across a broader interval.

`ps_density_plot()` estimates the positional distribution of the best hits for
a single motif.

```{r density-plot, fig.width=7.5, fig.height=4.6, dpi=120}
ps_density_plot(
    irf3,
    shift = -950, st = "strict",
    window = c(-950, 50), bins = 20
)
```

Higher-scoring IRF3 matches become more frequent towards the TSS, while a
substantial fraction remain distributed farther upstream.

The `window` argument specifies the promoter interval and enables reflection
boundary correction for the kernel density estimate. The boundaries correspond
to valid motif start positions, accounting for motif width. IRF3 is
`r ncol(irf3)` bp wide, so its curve stops `r ncol(irf3) - 1` bp short of the
downstream end.

The `bins` argument adds a histogram on the same density scale. Unlike the
kernel estimate, the histogram does not use a smoothing bandwidth, although
its appearance depends on bin width and placement. In this example, the final
bin is approximately 1.6 times as high as the smoothed curve, illustrating how
smoothing can reduce the apparent height of a narrow peak.

The `st` argument selects the displayed hits: `"all"` includes every promoter,
`"loose"` selects scores at or above the background mean, and `"strict"` selects
scores at least one background standard deviation above that mean. This can
highlight positional patterns among stronger sequence matches.

```{r density-counts, include=FALSE}
irf3_scores <- ps_hits_score(irf3)
loose_n <- sum(irf3_scores >= ps_bg_avg(irf3))
strict_n <- sum(irf3_scores >= ps_bg_avg(irf3) + ps_bg_std_dev(irf3))
```

For IRF3, the loose threshold retains `r loose_n` of the `r ps_hits_size(irf3)`
promoters and the strict threshold retains `r strict_n`. These display filters
are applied after scanning and do not alter the enrichment statistic, which
uses all scanned promoters.

`ps_hit_score_plot()` plots best-hit score against position for each promoter
and shows the `st` score thresholds as reference lines.

```{r hit-score-plot, fig.width=9.5, fig.height=4.6, dpi=115}
ps_hit_score_plot(
    results_subset[c("MA1418.1", "MA1513.1")],
    shift = -950
)
```

The joint display distinguishes positional concentration from match strength.
Many KLF15 hits in the final 200 bp approach the maximum score. IRF3 hits show
a weaker positional concentration and generally remain closer to the score
threshold. The combined score and position view therefore provides information
that a positional density alone does not convey.

# Analysing a motif-selected promoter subset

The preceding analysis considers all retina-only upregulated promoters. A
complementary analysis restricts the foreground to promoters containing a
high-scoring match to a selected motif.

`pscan_filtered()` scans the promoters with one chosen matrix, keeps those
scoring at least `n` standard deviations above that motif's background mean, and
rescans the retained promoters against the supplied collection.

The following code can be run after the preceding data-preparation steps. It is
not evaluated during vignette compilation to limit computation time.

```{r pscan-filtered, eval=FALSE}
irf3_matrix <- background[[which(TFBSTools::ID(background) == "MA1418.1")]]

filtered <- pscan_filtered(
    retina_up_promoters, irf3_matrix, n = 1,
    background = background,
    BPPARAM = BiocParallel::SerialParam()
)
ps_results_table(filtered)
```

Filtering by IRF3 score retains `r precomputed$filtered_irf3$n_retained` of the
`r precomputed$filtered_irf3$n_input` retina-only upregulated promoters.

```{r filtered-table, echo=FALSE}
ft <- precomputed$filtered_irf3$table
ft <- head(ft[order(-ft$ZSCORE), ], 12)
knitr::kable(
    data.frame(
        motif_id = ft$motif_id, NAME = ft$NAME,
        ZSCORE = round(ft$ZSCORE, 2),
        P.VALUE = format(ft$P.VALUE, scientific = TRUE, digits = 3)
    ),
    row.names = FALSE, align = c("l", "l", "r", "r"),
    caption = paste(
        "Enrichment within the", precomputed$filtered_irf3$n_retained,
        "promoters with a high-scoring IRF3 match."
    )
)
```

Interferon-related motifs have the highest ranks after selection by IRF3 score.
The GC-rich zinc-finger motifs that ranked highest in the unfiltered analysis
are absent from the leading results shown here.

```{r filtered-contrast, echo=FALSE}
zinc <- c("MA1513.1", "MA1650.1", "MA0162.4")
full <- precomputed$filtered_irf3$table
knitr::kable(
    data.frame(
        NAME = enrichment[zinc, "NAME"],
        `Z all promoters` = round(enrichment[zinc, "ZSCORE"], 2),
        `Z IRF3-positive` = round(
            full$ZSCORE[match(zinc, full$motif_id)], 2
        ),
        check.names = FALSE
    ),
    row.names = FALSE,
    caption = paste(
        "GC-rich motif scores before and after promoter selection",
        "by IRF3 score."
    )
)
```

This change is consistent with heterogeneity in the promoter set, as suggested
by the score heatmap and position map. These analyses use the same sequence
scores and are complementary descriptions, not independent validation.

Filtering provides a way to examine subsets of a heterogeneous foreground.
Because selection is based on a motif score, enrichment of that motif and
similar motifs is expected and should not be treated as independent evidence
of regulation.

# Summary

The core workflow comprises four steps: select a background, map genes to
promoters with a matching definition, scan the sequences, and inspect the
motif enrichment table.

```{r workflow, eval=FALSE}
background <- ps_retrieve_bg(
    "Jaspar2020", "mm", c(-950, 50), "mm10", version = 2
)
retina_up_promoters <- ps_select_promoters(
    genes = retina_up_genes, promoter_sequences = promoters_mm10,
    annotation = annotation, gene_col = "ENSEMBL", transcript_col = "REFSEQ",
    mode = "representative", return = "sequences"
)
results <- pscan(retina_up_promoters, background)
ps_results_table(results)
```

Additional analyses examine per-promoter scores (`ps_z_table()`), match
positions and sequences (`ps_hits_table()`, `ps_hitpos_map()`,
`ps_density_plot()`, `ps_hit_score_plot()`), and promoter subsets selected by
motif score (`pscan_filtered()`). These complementary summaries support the
interpretation of candidate regulatory patterns.

# Appendix: Comparison with the published analysis

The source study used the original Pscan implementation. The comparison above
follows the layout of Supplementary Figure S2C. In this PscanR analysis, all
nine IRF matrices in JASPAR 2020 meet an FDR threshold of 0.05 in the
retina-only upregulated set, together with `STAT1::STAT2`, `Stat2`, and `Stat4`.
None meets this threshold in the other five sets.

Two methodological differences affect the comparison.

```{r selection-counts, include=FALSE}
retina_selection <- ps_selection_summary(retina_mapping)$overall
```

**Transcript selection.** The published analysis does not document which
transcript was used for a gene with several. `ps_select_promoters()` makes that
choice explicit and deterministic, and records it: of the
`r retina_selection$genes_mapped` genes mapped here,
`r retina_selection$genes_with_choice` had more than one candidate. Exact motif
ranks may therefore differ from the original web-server run.

**Multiple testing.** PscanR reports nominal p-values and Benjamini–Hochberg
adjusted p-values in the `FDR` column; the published figure marks nominal
thresholds. Alternative corrections can be applied to the nominal p-values
in the complete result table, retaining the appropriate testing family.

# 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()
```