Motif enrichment analysis with PscanR: mouse retinal promoters
2026-11-11
Contents
- 1 Installation
- 2 Introduction
- 3 The background
- 4 From genes to promoters
- 5 A small executable scan
- 6 Running the scan
- 7 Checking the result against a second gene set
- 8 Inspecting the result object
- 9 Motif score patterns across promoters
- 10 Positional distribution of motif matches
- 11 Analysing a motif-selected promoter subset
- 12 Summary
- 13 Appendix: Comparison with the published analysis
- 14 Data provenance
- 15 Session information
1 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.
2 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.
2.1 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
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.
membership <- read.delim(
system.file(
"extdata", "vignettes", "mouse_hd_retina",
"paper_subset_membership.tsv",
package = "PscanR"
),
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| set | definition | genes |
|---|---|---|
| a | retina-only down | 935 |
| b | shared down | 143 |
| c | striatum-only down | 620 |
| d | retina-only up | 550 |
| e | shared up | 25 |
| f | striatum-only up | 151 |
The rest of this vignette follows set d, the
550 genes upregulated in retina and not in striatum, and
uses set f, the 151 genes upregulated in striatum and not
in retina, as the comparison.
3 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."
)| organism | assembly | upstream | downstream | jaspar_release | motif_count | promoter_count |
|---|---|---|---|---|---|---|
| mm | mm10 | 950 | 50 | 2020 | 746 | 31504 |
ps_retrieve_bg() retrieves the background statistics
and associates them with the matching JASPAR matrices in a PSMatrixList.
background <- ps_retrieve_bg(
jaspar_release = "Jaspar2020",
org = "mm",
prom_reg = c(-950, 50),
assembly = "mm10",
version = 2
)
length(background)
#> [1] 746This 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.
4 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(system.file(
"extdata", "vignettes", "mouse_hd_retina", "promoter_selection_example.rds",
package = "PscanR"
))
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 matchingThe returned mapping records the selected promoter and the selection rule:
| gene | promoter_id | scheme | suffix_role | transcript_class | selection_source | n_candidates | decided_by |
|---|---|---|---|---|---|---|---|
| ENSMUSG00000000131 | NM_001311143.2 | refseq | version | NM | RefSeq NM | 2 | tie_break |
| ENSMUSG00000000276 | NM_001362845.1 | refseq | version | NM | RefSeq NM | 3 | tie_break |
| ENSMUSG00000000305 | NM_001316723.1 | refseq | version | NM | RefSeq NM | 2 | tie_break |
| ENSMUSG00000000308 | NM_001355069.1 | refseq | version | NM | RefSeq NM | 2 | tie_break |
| ENSMUSG00000000378 | NM_001190343.1 | refseq | version | NM | RefSeq NM | 3 | tie_break |
| ENSMUSG00000000420 | NM_001160404.1 | refseq | version | NM | RefSeq NM | 3 | tie_break |
| ENSMUSG00000000530 | NM_001277255.1 | refseq | version | NM | RefSeq NM | 5 | tie_break |
| ENSMUSG00000000555 | NM_001314041.1 | refseq | version | NM | RefSeq NM | 2 | tie_break |
| ENSMUSG00000053420 | NR_033209.1 | refseq | version | NR | RefSeq NR | 1 | rule |
| ENSMUSG00000054057 | NR_028376.1 | refseq | version | NR | RefSeq NR | 2 | tie_break |
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(system.file(
"extdata", "vignettes", "mouse_hd_retina",
"mm10_950u_50d_refseq_snapshot.rds",
package = "PscanR"
))
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(system.file(
"extdata", "vignettes", "mouse_hd_retina", "uc2_vignette_data.rds",
package = "PscanR"
))
retina_mapping <- precomputed$retina_mappingps_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 2985Nearly 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 5365 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.45198376 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.1The 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-09GC-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.03615341BG_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")
7 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 22The 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.
7.1 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.
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 "" "" "" "" ...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.
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.
| set | nominal P < 0.05 | FDR < 0.05 | mean Z |
|---|---|---|---|
| a: retina-only down | 0 | 0 | -2.31 |
| b: shared down | 0 | 0 | -1.75 |
| c: striatum-only down | 0 | 0 | -4.10 |
| d: retina-only up | 14 | 12 | 3.22 |
| e: shared up | 3 | 0 | 0.75 |
| f: striatum-only up | 0 | 0 | -0.80 |
Set e contains only 24 promoters. Its small size limits
interpretation of the normal approximation
used for enrichment.
8 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] 536The 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 + AAGAAATAGAAACCGAAACTGThe 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"9 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.31ps_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.
10 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.
11 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.
| motif_id | NAME | ZSCORE | P.VALUE |
|---|---|---|---|
| MA1418.1 | IRF3 | 20.98 | 5.20e-98 |
| MA0653.1 | IRF9 | 16.35 | 2.17e-60 |
| MA1419.1 | IRF4 | 13.75 | 2.45e-43 |
| MA0652.1 | IRF8 | 13.59 | 2.16e-42 |
| MA0517.1 | STAT1::STAT2 | 12.75 | 1.52e-37 |
| MA0051.1 | IRF2 | 12.48 | 4.84e-36 |
| MA0050.2 | IRF1 | 11.36 | 3.13e-30 |
| MA0772.1 | IRF7 | 10.89 | 6.46e-28 |
| MA1420.1 | IRF5 | 10.78 | 2.05e-27 |
| MA1509.1 | IRF6 | 7.87 | 1.82e-15 |
| MA1623.1 | Stat2 | 7.79 | 3.30e-15 |
| MA0518.1 | Stat4 | 4.20 | 1.32e-05 |
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.
| NAME | Z all promoters | Z IRF3-positive |
|---|---|---|
| KLF15 | 8.27 | 2.51 |
| ZBTB14 | 7.95 | 1.79 |
| EGR1 | 6.67 | 1.81 |
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.
12 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.
13 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.
14 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.
15 Session information
sessionInfo()
#> R version 4.6.1 (2026-06-24)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 24.04.5 LTS
#>
#> Matrix products: default
#> BLAS: /home/biocbuild/bbs-3.24-bioc/R/lib/libRblas.so
#> LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.12.0 LAPACK version 3.12.0
#>
#> locale:
#> [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
#> [3] LC_TIME=en_GB LC_COLLATE=C
#> [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: America/New_York
#> tzcode source: system (glibc)
#>
#> attached base packages:
#> [1] stats graphics grDevices utils datasets methods base
#>
#> other attached packages:
#> [1] PscanRBackgrounds_0.99.7 PscanR_0.99.4 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 Biostrings_2.81.9
#> [7] S7_0.2.2 bitops_1.1-0
#> [9] fastmap_1.2.0 RCurl_1.98-1.20
#> [11] BiocFileCache_3.3.3 GenomicAlignments_1.49.2
#> [13] XML_3.99-0.25 digest_0.6.39
#> [15] DirichletMultinomial_1.55.1 lifecycle_1.0.5
#> [17] pwalign_1.9.1 KEGGREST_1.53.6
#> [19] RSQLite_3.53.3 magrittr_2.0.5
#> [21] compiler_4.6.1 rlang_1.3.0
#> [23] sass_0.4.10 tools_4.6.1
#> [25] yaml_2.3.12 rtracklayer_1.73.0
#> [27] knitr_1.52 labeling_0.4.3
#> [29] bit_4.6.0 curl_8.1.0
#> [31] RColorBrewer_1.1-3 BiocParallel_1.47.0
#> [33] withr_3.0.3 purrr_1.2.2
#> [35] BiocGenerics_0.59.13 grid_4.6.1
#> [37] stats4_4.6.1 ExperimentHub_3.3.2
#> [39] caTools_1.18.4 ggplot2_4.0.3
#> [41] scales_1.4.0 gtools_3.9.5
#> [43] tinytex_0.61 dichromat_2.0-1
#> [45] cli_3.6.6 rmarkdown_2.32
#> [47] crayon_1.5.3 generics_0.1.4
#> [49] otel_0.3.0 TFMPvalue_1.0.0
#> [51] httr_1.4.9 rjson_0.2.23
#> [53] BiocBaseUtils_1.15.1 DBI_1.3.0
#> [55] cachem_1.1.0 parallel_4.6.1
#> [57] AnnotationDbi_1.75.2 BiocManager_1.30.27
#> [59] XVector_0.53.0 restfulr_0.0.17
#> [61] matrixStats_1.5.0 vctrs_0.7.3
#> [63] jsonlite_2.0.0 bookdown_0.48
#> [65] IRanges_2.47.5 S4Vectors_0.51.10
#> [67] bit64_4.8.6 magick_2.9.1
#> [69] seqLogo_1.79.0 jquerylib_0.1.4
#> [71] glue_1.8.1 TFBSTools_1.51.0
#> [73] codetools_0.2-20 gtable_0.3.6
#> [75] BiocVersion_3.24.0 BiocIO_1.23.3
#> [77] GenomicRanges_1.65.4 tibble_3.3.1
#> [79] pillar_1.11.1 rappdirs_0.3.4
#> [81] htmltools_0.5.9 Seqinfo_1.3.2
#> [83] JASPAR2020_0.99.10 BSgenome_1.81.1
#> [85] httr2_1.3.0 dbplyr_2.6.0
#> [87] R6_2.6.1 evaluate_1.0.5
#> [89] Biobase_2.73.2 AnnotationHub_4.3.3
#> [91] png_0.1-9 Rsamtools_2.29.0
#> [93] cigarillo_1.3.1 pheatmap_1.0.13
#> [95] memoise_2.0.1 bslib_0.12.0
#> [97] Rcpp_1.1.2 xfun_0.62
#> [99] pkgconfig_2.0.3