---
title: "Motif enrichment analysis with PscanR: Arabidopsis BPC target promoters"
author:
    - Federico Zambelli
date: "`r Sys.Date()`"
output:
    BiocStyle::html_document:
        toc: true
        toc_depth: 2
vignette: >
    %\VignetteIndexEntry{4. PscanR: Arabidopsis BPC target 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", "arabidopsis_bpc", ...,
        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 *A. thaliana*
promoters and the functions available for examining the results.

## The data

The example uses gene expression data from a study of BASIC PENTACYSTEINE (BPC)
transcription factors, which control differentiation of the medial domain of
the *Arabidopsis* gynoecium.

> Caselli F, Palermiti M, Petrella R, Morlacchi VA, Dünser K, Kleine-Vehn J,
> Chiara M, Gregis V. *A novel role of BPCs in the control of medial domain
> differentiation during gynoecium development in Arabidopsis thaliana.*
> Plant Mol Biol 116(1):4 (2025).
> [doi:10.1007/s11103-025-01662-x](https://doi.org/10.1007/s11103-025-01662-x)

The input is a two-column table containing gene identifiers and directions of
expression change.

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

```{r read-gene-sets}
gene_sets <- read.delim(
    fixture("arabidopsis_bpc_gene_sets.tsv"),
    stringsAsFactors = FALSE
)
str(gene_sets)
table(gene_sets$direction)
```

The rest of this vignette follows the `r sum(gene_sets$direction == "down")`
genes with reduced expression when BPC function is impaired.

# The background

Precomputed backgrounds are indexed by organism, genome assembly, promoter
window, and motif collection. `ps_available_bg()` lists the available resources.

```{r discover-background}
ps_available_bg(keyword = "TAIR9", details = TRUE)
```

`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 = "at",
    prom_reg = c(-1000, 0),
    assembly = "TAIR9",
    version = 2,
    source = "zenodo"
)
length(background)
```

This returns 530 matrices from the JASPAR 2020 CORE plant 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 1000 bp upstream 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 identifiers or 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 following minimal example illustrates TAIR transcript selection:

```{r promoter-scheme, message=TRUE}
promoters <- c("AT1G01110.2", "AT1G01110.1", "AT1G01160.1")
annotation <- data.frame(
    GENE = sub("\\..*$", "", promoters),
    TX = promoters
)

selection <- ps_select_promoters(
    genes = c("AT1G01110", "AT1G01160"),
    promoter_ids = promoters,
    annotation = annotation,
    gene_col = "GENE", transcript_col = "TX",
    mode = "representative"
)
```

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"
    )],
    caption = "Selected promoters and the rules used for their selection."
)
```

For *Arabidopsis* the rule is to prefer the `.1` primary gene model and fall
back to the lowest available variant.

This dataset requires promoter sequences matching the background selected
above: TAIR9 annotation, 1000 bp upstream of the TSS.

Annotations derived from newer gene models should not be assumed to reproduce
TAIR9 promoters, even when transcript identifiers are shared. Changes to
annotated transcript ends can shift the TSS and the extracted sequence. This
example therefore uses the TAIR9 annotation specified by the background.

The following code illustrates promoter extraction and selection. It is not
evaluated during vignette compilation because it requires an external
annotation file and the `txdbmaker`, `GenomicFeatures`, `GenomicRanges`,
`GenomeInfoDb`, `Biostrings`, and `BSgenome.Athaliana.TAIR.TAIR9` packages.

```{r tair9-url, include=FALSE}
tair9_url <- paste0(
    "https://www.arabidopsis.org/download/list?",
    "dir=Genes/TAIR9_genome_release/tair9_gff3"
)
```

`TAIR9_GFF3_genes.gff` can be downloaded from the
[TAIR9 GFF3 release directory](`r tair9_url`).
Set `tair9_gff` to the local path of the downloaded file.

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

down_genes <- gene_sets$gene_id[gene_sets$direction == "down"]

# Use the TAIR9 annotation used to construct the background.

tair9_gff <- "TAIR9_GFF3_genes.gff" # path to the downloaded annotation

txdb <- txdbmaker::makeTxDbFromGFF(
    tair9_gff, format = "gff3",
    dataSource = "TAIR9", organism = "Arabidopsis thaliana"
)
genome <- BSgenome.Athaliana.TAIR.TAIR9::BSgenome.Athaliana.TAIR.TAIR9

# The background covers the five nuclear chromosomes.
tx <- GenomicFeatures::transcripts(txdb, columns = "tx_name")
tx <- tx[as.character(GenomicRanges::seqnames(tx)) %in% paste0("Chr", 1:5)]
names(tx) <- GenomicRanges::mcols(tx)$tx_name

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

# Exclude out-of-bounds ranges before sequence extraction; getSeq() otherwise
# raises an error. One transcript is affected here.
inside <- GenomicRanges::start(ranges) >= 1L &
    GenomicRanges::end(ranges) <= GenomeInfoDb::seqlengths(genome)[
        as.character(GenomicRanges::seqnames(ranges))
    ]
ranges <- ranges[inside]

promoters_tair9 <- Biostrings::getSeq(genome, ranges)
names(promoters_tair9) <- names(ranges)

# Exclude one additional promoter with more than 50% ambiguous N bases.
promoters_tair9 <- promoters_tair9[
    Biostrings::letterFrequency(promoters_tair9, "N", as.prob = TRUE) <= 0.5
]

annotation <- data.frame(
    GENE = sub("\\..*$", "", names(promoters_tair9)),
    TX = names(promoters_tair9)
)

down_mapping <- ps_select_promoters(
    genes = down_genes,
    promoter_sequences = promoters_tair9,
    annotation = annotation,
    gene_col = "GENE", transcript_col = "TX",
    mode = "representative",
    return = "mapping"
)

down_promoters <- promoters_tair9[down_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("arabidopsis_bpc_results.rds"))
down_mapping <- precomputed$down_mapping
```

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

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

Unmapped genes have no eligible promoter in the supplied reference. In this
dataset, annotation coverage accounts for several missing genes, as discussed
in the appendix.
`selection_source` records the rule used to select each promoter:

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

# A small executable scan

The following example uses twelve genes from the same downregulated list and
all their candidate TAIR9 promoters, extracted over -1000/0 from the archived
annotation snapshot. Selection preserves splice-variant suffixes. Two motifs
keep execution short; this demonstration does not replace the full analysis,
and its FDR adjustment covers only those two motifs.

```{r small-scan}
small_example <- readRDS(fixture("promoter_selection_example.rds"))
small_sequences <- ps_select_promoters(
    genes = small_example$genes,
    promoter_sequences = small_example$promoter_sequences,
    annotation = small_example$annotation, gene_col = "GENE",
    transcript_col = "TX", scheme = "tair", mode = "representative",
    return = "sequences"
)
small_background <- precomputed$results_subset[1:2]
small_results <- pscan(small_sequences, small_background)
ps_results_table(small_results)
```

# Running the scan

`pscan()` takes foreground promoter sequences and a background `PSMatrixList`,
and returns a `PSMatrixList` containing the per-promoter results for each motif.

The `BPPARAM` argument configures execution through BiocParallel:

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

- `SnowParam()` uses independent R workers connected through sockets. It
    supports Windows and Unix-like systems, with additional startup and
    data-transfer 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.

The following scan is not evaluated during vignette compilation to limit
computation time. It can be run after the preceding data-preparation steps;
runtime depends on the hardware and execution strategy.

```{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(
    down_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 downregulated promoters.
enrichment <- precomputed$enrichment_down
dim(enrichment)

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

The subset includes BPC motifs relevant to the original study, the
highest-ranked G-box-associated motifs, and two unrelated motifs 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), ], 6)
```

G-box-associated motifs from the bHLH and bZIP families occupy the highest
ranks. ABI5, ABF1, and ABF2 are bZIP factors whose motifs contain an ACGT core.

The background parameters remain associated with the result and can be retrieved
from it at any point with `ps_get_bg_table()`:

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

`BG_SIZE` is the number of *Arabidopsis* 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 the G-box binders rank above the three BPC motifs, as they do
in the full table above.

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

```{r barplot-composable, eval=FALSE}
ps_motif_barplot(results_subset, group = "class") +
    ggplot2::labs(title = "Downregulated promoters")
```

# 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; `MA0931.1` selects
the ABI5 entry:

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

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

`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(abi5), 5)
```

The ABI5 matches illustrate sequence variation around the G-box core `CACGTG`.
Oligonucleotides are reported in the forward orientation of the supplied
promoter sequence, including matches assigned to the reverse strand.

# 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 = 0.01,
    show_rownames = FALSE, treeheight_row = 0
)
```

The column dendrogram groups the three BPC motifs, consistent with their
similar GA-rich sequence preferences. A separate group contains
G-box-associated motifs, including HY5, ABF/AREB bZIPs, and PIF and SPT bHLHs.
The row-scaled heatmap suggests that promoters with relatively high BPC scores
differ from those with relatively high G-box motif scores. This heterogeneity
is not apparent from the motif-level enrichment table alone.

# Positional distribution of motif matches

`ps_hitpos_map()` displays the position of each promoter's best motif match.

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

`shift = -1000` relabels the axis so that coordinates are relative to the
transcription start site, matching how the promoters were extracted.

`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(
    abi5,
    shift = -1000, st = "strict",
    window = c(-1000, 0), bins = 20
)
```

Higher-scoring ABI5 matches are concentrated near the TSS, with the estimated
density peaking within approximately 100 bp 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.

The `bins` argument adds a histogram of the same positions on the density
scale, allowing comparison with the smoothed curve.

The BPC6 motif, with its GA-rich sequence preference, provides a comparison
with the ACGT-containing ABI5 motif.

```{r density-plot-bpc6, fig.width=7.5, fig.height=4.6, dpi=120}
bpc6_hits <- results_subset[["MA1402.1"]]

ps_density_plot(
    bpc6_hits,
    shift = -1000, st = "strict",
    window = c(-1000, 0), bins = 20
)
```

Both profiles show concentration near the TSS. The ABI5 density peaks
approximately 75–100 bp upstream and then decreases, whereas the BPC6 profile
increases towards the final bin. BPC6 has width `r ncol(bpc6_hits)` bp, so the
last valid motif start lies `r ncol(bpc6_hits) - 1` bp before the end of the
window.

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 matches. These display filters
do not alter the motif enrichment statistics.

`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("MA1402.1", "MA0931.1")],
    shift = -1000
)
```

Higher-scoring BPC6 hits are concentrated within the final few hundred bases,
whereas lower-scoring hits are distributed across the window. Higher-scoring
ABI5 hits show a broader distribution and generally occur farther upstream.
The joint display distinguishes positional concentration from match strength.

# Analysing a motif-selected promoter subset

The preceding analysis considers all downregulated 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}
bpc6 <- background[[which(TFBSTools::ID(background) == "MA1402.1")]]

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

Filtering by BPC6 score retains `r precomputed$filtered_bpc6$n_retained` of the
`r precomputed$filtered_bpc6$n_input` downregulated promoters.

```{r filtered-table, echo=FALSE}
ft <- precomputed$filtered_bpc6$table
ft <- head(ft[order(-ft$ZSCORE), ], 10)
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_bpc6$n_retained,
        "promoters with a high-scoring BPC6 match."
    )
)
```

The GA-rich BPC and RAMOSA1 motifs have the highest ranks after selection by
BPC6 score. G-box-associated bZIP and bHLH motifs are absent from the leading
results shown here, which instead include C2H2 and Dof motifs. This change is
consistent with different sequence properties in the BPC6-selected subset.

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", "at", c(-1000, 0), "TAIR9", version = 2
)
down_promoters <- ps_select_promoters(
    genes = down_genes, promoter_sequences = promoters_tair9,
    annotation = annotation, gene_col = "GENE", transcript_col = "TX",
    mode = "representative", return = "sequences"
)
results <- pscan(down_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 reports Pscan results that permit motif-level comparison.
Using the comparison criterion described below, this analysis recovers 164 of
the 168 motifs listed for the published downregulated set. The agreement is
substantial but incomplete. The following methodological details affect the
comparison.

**The promoter window.** The published analysis used −1000/+100 relative to the
TSS. The bundled catalogue does not include that window, so this analysis uses
−1000/0. A separate analysis of the same genes with the −950/+50 background
gave a Spearman correlation of ρ = 0.99 across all 530 motif z-scores, although
some significance classifications changed. This comparison supports similarity
between the two tested windows; it does not directly measure the effect of
adding the full 100 bp downstream interval used in the publication.

**The annotation release.** The published promoters were defined on TAIR10 gene
models; this background uses TAIR9. Shared genome coordinates do not guarantee
identical transcript models or TSS annotations. In the recorded reference
comparison, 2244 of the study's 2253 genes have promoters in the TAIR9 universe
(99.6%). The nine missing genes include loci absent from the older annotation
and a mitochondrial gene outside the nuclear reference set; they occur in both
expression directions. High gene coverage does not establish promoter-level
identity between annotation releases.

**Multiple testing correction.** PscanR reports nominal p-values and
Benjamini–Hochberg adjusted p-values in the `FDR` column. The
[published methods](https://doi.org/10.1007/s11103-025-01662-x) specify
Benjamini–Hochberg correction with a 0.01 threshold. However, inspection of the
deposited Supplementary Table 3 in the local reproduction found adjusted values
approximately equal to nominal p-values multiplied by 530, with listed values
below 0.1.
The overlap of 164 out of 168 therefore refers to the reproduction's
Bonferroni criterion, `p * 530 < 0.1`, chosen to match those deposited values;
it is not the count obtained from the default PscanR `FDR` column. 

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