Installation

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

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.

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

membership <- read.delim(
    fixture("paper_subset_membership.tsv"),
    stringsAsFactors = FALSE
)
str(membership)
#> 'data.frame':    2424 obs. of  2 variables:
#>  $ subset         : chr  "a" "a" "a" "a" ...
#>  $ ensembl_gene_id: chr  "ENSMUSG00000060890" "ENSMUSG00000058831" "ENSMUSG00000030324" "ENSMUSG00000009108" ...
gene_sets <- split(membership$ensembl_gene_id, membership$subset)
lengths(gene_sets)
#>   a   b   c   d   e   f 
#> 935 143 620 550  25 151

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:

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.

background <- ps_retrieve_bg(
    JASPAR_matrix = "Jaspar2020",
    org = "mm",
    prom_reg = c(-950, 50),
    assembly = "mm10",
    version = 2,
    source = "zenodo"
)
length(background)
#> [1] 746

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:

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"
)
#> Detected transcript scheme: "refseq" (25/25 identifiers matched)
#>   NM_001355069.1 -> transcript NM_001355069
#>   the trailing .N is a VERSION and is normalised away for matching

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

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.

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.

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.

ps_selection_summary(retina_mapping)$overall
#>   scheme genes_requested genes_mapped genes_unmapped genes_with_choice
#> 1 refseq             550          536             14               260
#>   genes_decided_by_tie_break candidates_excluded_protein
#> 1                        248                        2985

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:

head(attr(retina_mapping, "ps_unmapped_genes"), 4)
#> [1] "ENSMUSG00000086583" "ENSMUSG00000064357" "ENSMUSG00000058443"
#> [4] "ENSMUSG00000064356"

selection_source records the rule used to select each promoter:

ps_selection_summary(retina_mapping)$by_source
#>   selection_source promoters
#> 1        RefSeq NM       536

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.

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)
#>          NAME    BG_AVG   BG_STDEV    FG_AVG    ZSCORE   P.VALUE       FDR
#> MA0652.1 IRF8 0.7459045 0.04172396 0.7482274 0.1760553 0.4301252 0.4519837
#> MA1418.1 IRF3 0.7357228 0.03982087 0.7372421 0.1206510 0.4519837 0.4519837

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.

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

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.

enrichment <- ps_results_table(results)

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

# Saved enrichment table for retina-only upregulated promoters.
enrichment <- precomputed$enrichment_d
dim(enrichment)
#> [1] 746   7

# A motif subset of the PSMatrixList retains the per-promoter results.
results_subset <- precomputed$results_subset_d
results_subset
#> PSMatrixList of length 8
#> names(8): MA1418.1 MA0652.1 MA0517.1 MA1623.1 MA1513.1 MA1650.1 MA0162.4 MA0668.1

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.

TFBSTools::name(results_subset)
#>       MA1418.1       MA0652.1       MA0517.1       MA1623.1       MA1513.1 
#>         "IRF3"         "IRF8" "STAT1::STAT2"        "Stat2"        "KLF15" 
#>       MA1650.1       MA0162.4       MA0668.1 
#>       "ZBTB14"         "EGR1"      "NEUROD2"

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

head(enrichment[order(-enrichment$ZSCORE), ], 8)
#>            NAME    BG_AVG   BG_STDEV    FG_AVG   ZSCORE      P.VALUE
#> MA1513.1  KLF15 0.8662098 0.08924349 0.8981071 8.274845 6.431123e-17
#> MA1650.1 ZBTB14 0.8220133 0.06587730 0.8446258 7.946832 9.567053e-16
#> MA1564.1    SP9 0.8777526 0.05873299 0.8962573 7.294299 1.501093e-13
#> MA1584.1   ZIC5 0.7796885 0.04413328 0.7931302 7.051317 8.861623e-13
#> MA0732.1   EGR3 0.7905038 0.05461574 0.8066153 6.829679 4.255256e-12
#> MA0746.2    SP3 0.8646196 0.05901386 0.8819176 6.786167 5.757610e-12
#> MA0162.4   EGR1 0.8246556 0.06088805 0.8421905 6.667352 1.302301e-11
#> MA0697.1   ZIC3 0.7917738 0.04244414 0.8038072 6.563770 2.623201e-11
#>                   FDR
#> MA1513.1 4.797617e-14
#> MA1650.1 3.568511e-13
#> MA1564.1 3.732718e-11
#> MA1584.1 1.652693e-10
#> MA0732.1 6.348842e-10
#> MA0746.2 7.158628e-10
#> MA0162.4 1.387881e-09
#> MA0697.1 2.446135e-09

GC-rich motifs, including KLF15, ZBTB14, and SP and EGR zinc fingers, occupy the highest ranks. IRF3, the highest-ranked IRF motif, has rank 27 with a z-score of 5.44.

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

ps_get_bg_table(results_subset)
#>          BG_SIZE   BG_MEAN   BG_STDEV
#> MA1418.1   31504 0.7357228 0.03982087
#> MA0652.1   31504 0.7459045 0.04172396
#> MA0517.1   31504 0.8215885 0.04244822
#> MA1623.1   31504 0.8408313 0.04997154
#> MA1513.1   31504 0.8662098 0.08924349
#> MA1650.1   31504 0.8220133 0.06587730
#> MA0162.4   31504 0.8246556 0.06088805
#> MA0668.1   31504 0.8801757 0.03615341

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().

table(ps_motif_class(results_subset))
#> 
#> Basic helix-loop-helix factors (bHLH)              C2H2 zinc finger factors 
#>                                     1                                     3 
#>                   STAT domain factors            Tryptophan cluster factors 
#>                                     2                                     2

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 -4.2 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:

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.

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))
)
#>   set       definition promoters nominal fdr
#> 1   d   retina-only up       536     139 102
#> 2   f striatum-only up       146      72  22

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

ps_motif_barplot(precomputed$results_subset_f, group = "class")

The GC-rich zinc-finger motifs remain enriched: the KLF15 z-score changes from 8.27 to 4.01. In contrast, the IRF3 z-score changes from 5.44 to -1.40, 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.

irf_stat <- precomputed$irf_stat
str(irf_stat[c("subset", "motif", "ZSCORE", "significance")])
#> 'data.frame':    108 obs. of  4 variables:
#>  $ subset      : Factor w/ 6 levels "a: retina-only down",..: 1 1 1 1 1 1 1 1 1 1 ...
#>  $ motif       : Factor w/ 18 levels "Stat6 [MA0520.1]",..: 14 3 13 15 11 9 5 10 6 4 ...
#>  $ ZSCORE      : num  -0.678 -0.921 -1.076 -1.178 -1.354 ...
#>  $ significance: chr  "" "" "" "" ...

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"
    )

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:

irf3 <- results_subset[["MA1418.1"]]
TFBSTools::name(irf3)
#> [1] "IRF3"

ps_zscore(irf3)
#>        z 
#> 5.438571
ps_pvalue(irf3)
#> [1] 2.685481e-08
c(background = ps_bg_avg(irf3), foreground = ps_fg_avg(irf3))
#> background foreground 
#>  0.7357228  0.7450771
ps_hits_size(irf3)
#> [1] 536

The modest difference between foreground and background means produces an enrichment z-score of 5.44. This statistic concerns the mean across 536 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.

head(ps_hits_table(irf3, pos_shift = -950L), 5)
#>                                       SCORE  POS STRAND                 OLIGO
#> ENSMUSG00000040253|NM_001083312.2 0.9290642  -56      - CAGTTTCAGTTTCAGTTCTCC
#> ENSMUSG00000037921|NM_001293783.1 0.9244447  -61      + GATGAATGGAAACCGAAACTT
#> ENSMUSG00000034780|NM_001379324.1 0.9229804 -817      - AGGTTTCTGTTTCCGTTTCCT
#> ENSMUSG00000016921|NM_026499.4    0.9184677 -816      - TGCTTTCGGTTTCCTTACCTC
#> ENSMUSG00000046879|NM_001355757.1 0.9156476  -32      + AAGAAATAGAAACCGAAACTG

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:

pfm <- TFBSTools::Matrix(irf3)
paste(rownames(pfm)[apply(pfm, 2, which.max)], collapse = "")
#> [1] "GGGAAACGGAAACCGAAACTG"

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.

z <- ps_z_table(results_subset)
dim(z)
#> [1] 536   8
round(z[1:4, 1:5], 2)
#>                                MA1418.1 MA0652.1 MA0517.1 MA1623.1 MA1513.1
#> ENSMUSG00000060044|NM_177794.3    -0.37     0.08     1.36     1.32     1.26
#> ENSMUSG00000036030|NM_175485.4    -1.94    -0.39    -1.27    -1.45     1.31
#> ENSMUSG00000030111|NM_175628.3    -1.07    -1.42    -0.97    -0.81    -1.40
#> ENSMUSG00000032554|NM_133977.2     0.41     1.09     0.98    -0.63     1.31

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.

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.

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.

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 21 bp wide, so its curve stops 20 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.

For IRF3, the loose threshold retains 272 of the 536 promoters and the strict threshold retains 110. 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.

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.

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 110 of the 536 retina-only upregulated promoters.

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.

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.

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.

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 536 genes mapped here, 260 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

sessionInfo()
#> R version 4.6.1 (2026-06-24)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 26.04.1 LTS
#> 
#> Matrix products: default
#> BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
#> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.32.so;  LAPACK version 3.12.0
#> 
#> locale:
#>  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
#>  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
#>  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
#>  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
#>  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
#> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
#> 
#> time zone: Etc/UTC
#> tzcode source: system (glibc)
#> 
#> attached base packages:
#> [1] stats4    stats     graphics  grDevices utils     datasets  methods  
#> [8] base     
#> 
#> other attached packages:
#>  [1] org.Hs.eg.db_3.23.1  AnnotationDbi_1.75.2 Biobase_2.73.2      
#>  [4] Biostrings_2.81.9    Seqinfo_1.3.2        XVector_0.53.0      
#>  [7] IRanges_2.47.5       S4Vectors_0.51.10    BiocGenerics_0.59.12
#> [10] generics_0.1.4       PscanR_0.99.0        BiocStyle_2.41.0    
#> 
#> loaded via a namespace (and not attached):
#>  [1] tidyselect_1.2.1            dplyr_1.2.1                
#>  [3] farver_2.1.2                blob_1.3.0                 
#>  [5] filelock_1.0.3              S7_0.2.2                   
#>  [7] bitops_1.1-0                fastmap_1.2.0              
#>  [9] RCurl_1.98-1.20             BiocFileCache_3.3.0        
#> [11] GenomicAlignments_1.49.2    XML_3.99-0.25              
#> [13] digest_0.6.39               DirichletMultinomial_1.55.1
#> [15] lifecycle_1.0.5             pwalign_1.9.1              
#> [17] KEGGREST_1.53.6             RSQLite_3.53.3             
#> [19] magrittr_2.0.5              compiler_4.6.1             
#> [21] rlang_1.3.0                 sass_0.4.10                
#> [23] tools_4.6.1                 yaml_2.3.12                
#> [25] rtracklayer_1.73.0          knitr_1.52                 
#> [27] labeling_0.4.3              bit_4.6.0                  
#> [29] curl_8.0.0                  RColorBrewer_1.1-3         
#> [31] BiocParallel_1.47.0         withr_3.0.3                
#> [33] purrr_1.2.2                 sys_3.4.3                  
#> [35] grid_4.6.1                  caTools_1.18.4             
#> [37] ggplot2_4.0.3               scales_1.4.0               
#> [39] gtools_3.9.5                cli_3.6.6                  
#> [41] rmarkdown_2.32              crayon_1.5.3               
#> [43] otel_0.2.0                  TFMPvalue_1.0.0            
#> [45] httr_1.4.9                  rjson_0.2.23               
#> [47] BiocBaseUtils_1.15.1        DBI_1.3.0                  
#> [49] cachem_1.1.0                parallel_4.6.1             
#> [51] BiocManager_1.30.27         restfulr_0.0.17            
#> [53] matrixStats_1.5.0           vctrs_0.7.3                
#> [55] Matrix_1.7-6                jsonlite_2.0.0             
#> [57] bit64_4.8.6                 seqLogo_1.79.0             
#> [59] maketools_1.3.2             jquerylib_0.1.4            
#> [61] glue_1.8.1                  TFBSTools_1.51.0           
#> [63] codetools_0.2-20            gtable_0.3.6               
#> [65] BiocIO_1.23.3               GenomicRanges_1.65.4       
#> [67] tibble_3.3.1                pillar_1.11.1              
#> [69] htmltools_0.5.9             JASPAR2020_0.99.10         
#> [71] BSgenome_1.81.1             httr2_1.3.0                
#> [73] dbplyr_2.6.0                R6_2.6.1                   
#> [75] evaluate_1.0.5              lattice_0.23-1             
#> [77] png_0.1-9                   Rsamtools_2.29.0           
#> [79] cigarillo_1.3.1             pheatmap_1.0.13            
#> [81] memoise_2.0.1               bslib_0.12.0               
#> [83] Rcpp_1.1.2                  xfun_0.61                  
#> [85] buildtools_1.0.0            pkgconfig_2.0.3