PscanR is being prepared for Bioconductor submission. After acceptance, install it with the following command (not executed while building this vignette):
For pre-submission installation from source, see the package README.
PscanR provides access to precomputed backgrounds for five model organisms: human, mouse, and fruit fly (RefSeq annotation), Arabidopsis (TAIR annotation), and yeast (SGD annotation).
ps_available_bg() lists the available resources:
head(ps_available_bg(), 4)
#> [1] "J2020_dm6_200u_50d_UCSC.psbg2.txt" "J2020_dm6_450u_50d_UCSC.psbg2.txt"
#> [3] "J2020_dm6_500u_0d_UCSC.psbg2.txt" "J2020_dm6_950u_50d_UCSC.psbg2.txt"This vignette describes custom background construction. A custom background is required when the catalogue does not include the desired:
The example reconstructs an existing catalogue entry for Saccharomyces cerevisiae, using promoters from −950 to +50 bp relative to the transcription start site (TSS). This permits comparison with the published background.
The bundled timing and verification results were generated using the developers’ computational environment. The code below documents the construction procedure; expensive or externally dependent steps are not evaluated during vignette compilation.
The yeast example permits a relatively short demonstration. In the recorded run, background construction took 155 seconds on one core, for 6,571 promoters and 183 matrices. Runtime for other organisms depends on the number and length of promoters, the motif collection, and the hardware.
The following code builds both background representations from a subset of real yeast promoters and two JASPAR motifs. This restricted background is an API demonstration, not a replacement for the complete promoter universe used for enrichment. Both computations run during vignette builds.
example_sequences <- readRDS(fixture("rp_promoters.rds"))
n_example <- min(40L, length(example_sequences))
example_sequences <- example_sequences[seq_len(n_example)]
example_motifs <- JASPAR2020::JASPAR2020
example_motifs <- TFBSTools::getMatrixSet(
example_motifs, list(collection = "CORE", tax_group = "fungi")
)[1:2]
short_example <- ps_build_bg(example_sequences, example_motifs)
short_file <- tempfile(fileext = ".txt")
ps_write_bg_to_file(short_example, short_file)
short_restored <- ps_retrieve_bg_from_file(short_file, example_motifs)
stopifnot(isTRUE(all.equal(ps_get_bg_table(short_example),
ps_get_bg_table(short_restored))))
full_example <- ps_build_bg(example_sequences, example_motifs, fullBG = TRUE)
full_file <- tempfile(fileext = ".rds")
saveRDS(full_example, full_file, compress = "xz")
full_restored <- readRDS(full_file)
stopifnot(identical(full_example, full_restored))
query <- names(example_sequences)[1:4]
queried <- pscan_full_bg(query, full_restored)
ps_results_table(queried)
#> NAME BG_AVG BG_STDEV FG_AVG ZSCORE P.VALUE FDR
#> MA0266.1 ABF2 0.8852288 0.06052463 0.9241023 1.284553 0.09947426 0.1989485
#> MA0265.1 ABF1 0.8392502 0.03361023 0.8192210 -1.191850 0.88333990 0.8833399
unlink(c(short_file, full_file))A background is constructed from reference sequences and motif matrices. Here, the reference comprises annotated nuclear promoters, defined using the same window as the intended foreground: −950 to +50 bp relative to the TSS.
The reference is specified by three components:
sgdGene table,
provided by the TxDb.Scerevisiae.UCSC.sacCer3.sgdGene
package, whose version records the annotation release used;BSgenome.Scerevisiae.UCSC.sacCer3, whose chromosome names
must agree with the annotation’s;Changes to the assembly or annotation can alter promoter sequences and background statistics. Record the reference versions and use a consistent promoter definition for foreground and background.
# UCSC's sgdGene table, packaged and versioned.
library(TxDb.Scerevisiae.UCSC.sacCer3.sgdGene)
txdb <- TxDb.Scerevisiae.UCSC.sacCer3.sgdGene
# The background covers the 16 nuclear chromosomes.
tx <- GenomicFeatures::transcripts(txdb, columns = c("tx_id", "tx_name"))
nuclear <- paste0("chr", as.roman(seq_len(16)))
tx <- tx[as.character(GenomicRanges::seqnames(tx)) %in% nuclear]
names(tx) <- GenomicRanges::mcols(tx)$tx_name
# promoters() is strand-aware and anchored on the transcription start site.
ranges <- GenomicRanges::promoters(tx, upstream = 950, downstream = 50)
genome <- BSgenome.Scerevisiae.UCSC.sacCer3::BSgenome.Scerevisiae.UCSC.sacCer3
# Exclude out-of-bounds ranges before sequence extraction; getSeq() otherwise
# raises an error.
lengths <- GenomeInfoDb::seqlengths(genome)[
as.character(GenomicRanges::seqnames(ranges))
]
inside <- GenomicRanges::start(ranges) >= 1L &
GenomicRanges::end(ranges) <= lengths
ranges <- ranges[inside]
prom_seq <- Biostrings::getSeq(genome, ranges)
names(prom_seq) <- names(ranges)This produces 6,639 promoter sequences, one per retained transcript.
Retain the original sequence names and duplicate sequences when
passing the extracted promoters to ps_build_bg(). The
function deduplicates sequences, retains those of the maximum supplied
length, and removes those with more than 50% ambiguous N
bases. Supply the intended full-length promoter window; the function
does not infer it from genomic coordinates. Preserving names before
deduplication is particularly important for full backgrounds, where the
identifier mapping must retain equivalent transcripts. Out-of-bounds
genomic ranges must still be excluded before sequence extraction, as
shown above.
The background size is the number of distinct promoter sequences remaining after these operations: 6,571 in this example.
Supply the motif collection as a PFMatrixList, as
defined by TFBSTools. The catalogued backgrounds use JASPAR CORE
matrices restricted to the organism’s taxonomic group.
ps_build_bg() retains the best-scoring window for each
promoter–motif pair, considering both strands. It then computes the
background mean and standard deviation of these scores and records the
number of retained promoters.
| Step | One core | Eight cores | Speedup |
|---|---|---|---|
| ps_build_bg(), 183 matrices | 154.9 | 28.1 | 5.5 |
Using eight cores (BPPARAM set to
BiocParallel::MulticoreParam(workers = 8)) produced a
5.5-fold speedup for the full motif collection in the recorded run.
Parallelism is distributed over motifs. Actual speedup also depends on
promoter count and length, worker overhead, and available resources.
The following output shows the header, the first three motif records, and the file size in bytes:
bg_file <- fixture("J2020_sacCer3_950u_50d_UCSC.rebuild.txt")
cat(readLines(bg_file, n = 4), sep = "\n")
#> [SHORT TFBS MATRIX]
#> MA0265.1 6571 0.8608245107 0.0489984005788794
#> MA0266.1 6571 0.886681178483697 0.0605950135938033
#> MA0267.1 6571 0.920478384662466 0.0572881109107041
file.size(bg_file)
#> [1] 9323The text format contains a header followed by one row per matrix, with four fields: matrix identifier, promoter count, mean score, and score standard deviation. Storing summary statistics keeps standard background files small.
ps_retrieve_bg_from_file() associates these statistics
with the supplied matrices by identifier and returns a
PSMatrixList for use with pscan().
background <- ps_retrieve_bg_from_file(bg_file, motifs)
background
#> PSMatrixList of length 183
#> names(183): MA0265.1 MA0266.1 MA0267.1 MA0268.1 ... MA1435.1 MA1436.1 MA1437.1
head(ps_get_bg_table(background), 3)
#> BG_SIZE BG_MEAN BG_STDEV
#> MA0265.1 6571 0.8608245 0.04899840
#> MA0266.1 6571 0.8866812 0.06059501
#> MA0267.1 6571 0.9204784 0.05728811ps_get_bg_table() retrieves the three statistical
columns from the reloaded object. Text serialisation preserves their
values to the precision of the file format; it does not guarantee
bit-for-bit equality with the original in-memory floating-point
values.
The reconstructed and catalogued backgrounds use the same annotation, genome, promoter window, and matrices. The following table reports their agreement in the recorded verification run.
knitr::kable(
data.frame(
Check = c(
"Matrices compared",
"BG_SIZE equal",
"BG_MEAN agreeing to 10 decimal places",
"BG_STDEV agreeing to 10 decimal places"
),
Result = c(
verification$motif_count,
verification$size_equal,
sprintf("%d / %d", verification$mean_agree_10dp,
verification$motif_count),
sprintf("%d / %d", verification$stdev_agree_10dp,
verification$motif_count)
)
),
align = c("l", "r")
)| Check | Result |
|---|---|
| Matrices compared | 183 |
| BG_SIZE equal | TRUE |
| BG_MEAN agreeing to 10 decimal places | 183 / 183 |
| BG_STDEV agreeing to 10 decimal places | 183 / 183 |
All matrices agree in background size, and both score statistics agree to ten decimal places in this comparison.
Floating-point arithmetic and text serialisation can introduce small numerical differences. Agreement to ten decimal places is the criterion used for this verification, not a guarantee of exact byte-equality across platforms. Those tiny differences across platforms are usually negligible and do not impact PscanR results in any meaningful way.
A standard background stores three summary statistics per matrix. Foreground analysis therefore requires promoter sequences extracted using the matching reference and window.
PscanR also supports full backgrounds, built with
fullBG = TRUE. These retain the per-promoter scan results,
allowing subsequent foreground analyses to retrieve existing hits by
identifier. This avoids repeated DNA scanning and is useful when
analysing many gene sets against the same background.
A full background retains the best-hit score, position, strand, and matched oligonucleotide for every retained promoter and matrix. Its identifier legend maps original sequence names to the representatives retained after deduplication.
pscan_full_bg() accepts transcript identifiers matching
this legend, retrieves the stored hits, and calculates foreground
enrichment without rescanning DNA.
# The bundled fullbg_motifs.tsv records the fifteen selected matrices.
subset_ids <- read.delim(fixture("fullbg_motifs.tsv"))$MATRIX
full_bg <- ps_build_bg(
prom_seq, motifs[subset_ids],
BPPARAM = BiocParallel::SerialParam(), fullBG = TRUE
)Pass sequences before deduplication. The legend is constructed from the input names. Deduplicating beforehand prevents the function from retaining identifiers associated with the removed duplicate sequences.
Storage requirements. Full backgrounds can use sizeable storage up to hundreds of megabytes depending on the number of promoters and matrices. The bundled yeast demonstration contains 15 matrices and occupies approximately 692 KiB. Linear extrapolation by matrix count gives approximately 8 MiB for all 183 yeast matrices. Complete full-background collections are not maintaned nor bundled with the package.
Save a full background as an .rds file to preserve its
PSMatrixList, hits, and identifier legend.
The reported file sizes use compress = "xz". This
compression method can reduce storage requirements at the cost of
additional compression time.
Reloading restores the hits and legend used in the following foreground query.
full_bg <- readRDS(fixture("sacCer3_950u_50d_fullBG_15motifs.rds"))
TFBSTools::name(full_bg)
#> MA0359.1 MA0378.1 MA0295.1 MA0265.1 MA0363.1 MA0303.1 MA0331.1 MA0341.1
#> "RAP1" "SFP1" "FHL1" "ABF1" "REB1" "GCN4" "MCM1" "MSN2"
#> MA0342.1 MA0319.1 MA0299.1 MA0357.1 MA0393.1 MA0412.1 MA0324.1
#> "MSN4" "HSF1" "GAL4" "PHO4" "STE12" "UME6" "LEU3"The foreground example uses yeast cytoplasmic ribosomal protein genes to illustrate retrieval of a promoter set expected to contain Rap1 motif matches.
The set contains 132 cytoplasmic ribosomal protein genes selected
from org.Sc.sgd.db by standard names matching
^RP[LS][0-9] and mapped to SGD systematic identifiers. The
background contains 132 of these identifiers.
rp_promoters <- readRDS(fixture("rp_promoters.rds"))
rp_ids <- names(rp_promoters)
length(rp_ids)
#> [1] 132by_legend <- pscan_full_bg(rp_ids, full_bg)
#> Detected transcript scheme: "sgd" (6639/6639 identifiers matched)
#> YAL067W-A -> transcript YAL067W-A
#> these identifiers carry no version or variant suffixThe message reports the identifier scheme detected in the background.
enrichment <- ps_results_table(by_legend)
head(enrichment[, c("NAME", "ZSCORE", "P.VALUE", "FDR")], 5)
#> NAME ZSCORE P.VALUE FDR
#> MA0359.1 RAP1 20.865051 5.563764e-97 8.345646e-96
#> MA0378.1 SFP1 3.267663 5.421973e-04 4.066480e-03
#> MA0324.1 LEU3 2.602672 4.625022e-03 2.312511e-02
#> MA0299.1 GAL4 2.378410 8.693741e-03 3.260153e-02
#> MA0295.1 FHL1 1.499070 6.692777e-02 2.007833e-01RAP1 is the highest-ranked motif, with a z-score of 20.9. SFP1 ranks second in this restricted motif collection. The displayed table reports both nominal and adjusted p-values.
This workflow performs identifier lookup and retrieval of stored hits for each query, followed by enrichment calculation. It avoids promoter extraction and DNA rescanning. Query identifiers must match the background legend; gene symbols cannot be substituted unless they are the identifiers used by that background. Here, the SGD systematic identifiers name the stored promoters.
Record the following choices to support reproducibility:
TxDb if
Bioconductor ships one for the organism. Otherwise a UCSC table via
txdbmaker::makeTxDbFromUCSC(), or a GFF/GTF via
txdbmaker::makeTxDbFromGFF(); with those, record which
release was used and when it was retrieved.BSgenome package if
one exists; otherwise a FASTA read with
Biostrings::readDNAStringSet(). Record the assembly and
ensure that chromosome names agree with the annotation.scheme argument of
pscan_full_bg() supports "refseq",
"sgd", "tair", and "generic";
"auto" detects recognised conventions. Generic identifiers
are matched without suffix removal.
# 1. Extract promoters using the intended foreground window.
genome <- BSgenome.Scerevisiae.UCSC.sacCer3::BSgenome.Scerevisiae.UCSC.sacCer3
library(TxDb.Scerevisiae.UCSC.sacCer3.sgdGene)
txdb <- TxDb.Scerevisiae.UCSC.sacCer3.sgdGene
tx <- GenomicFeatures::transcripts(txdb, columns = c("tx_id", "tx_name"))
nuclear <- paste0("chr", as.roman(seq_len(16)))
tx <- tx[as.character(GenomicRanges::seqnames(tx)) %in% nuclear]
names(tx) <- GenomicRanges::mcols(tx)$tx_name
ranges <- GenomicRanges::promoters(tx, upstream = 950, downstream = 50)
lengths <- GenomeInfoDb::seqlengths(genome)[
as.character(GenomicRanges::seqnames(ranges))
]
inside <- GenomicRanges::start(ranges) >= 1L &
GenomicRanges::end(ranges) <= lengths
ranges <- ranges[inside]
prom_seq <- Biostrings::getSeq(genome, ranges)
# 2. Select motif matrices for the organism's taxonomic group.
motifs <- TFBSTools::getMatrixSet(
JASPAR2020::JASPAR2020, list(collection = "CORE", tax_group = "fungi")
)
# 3. Construct the standard background.
background <- ps_build_bg(
prom_seq, motifs, BPPARAM = BiocParallel::SerialParam()
)
# 4. Save the background statistics.
ps_write_bg_to_file(background, "my_background.txt")
background <- ps_retrieve_bg_from_file("my_background.txt", motifs)
# Optionally, retain all hits in a full background for identifier-based queries.
full_bg <- ps_build_bg(prom_seq, motifs, fullBG = TRUE)
saveRDS(full_bg, "my_full_background.rds", compress = "xz")
full_bg <- readRDS("my_full_background.rds")
# Query identifiers must match the background's promoter naming scheme.
# This example uses SGD systematic identifiers for cytoplasmic ribosomal
# protein genes. An external list can be read from a file with readLines().
gene_names <- unlist(as.list(org.Sc.sgd.db::org.Sc.sgdGENENAME))
my_gene_ids <- names(gene_names)[grepl("^RP[LS][0-9]", gene_names)]
ps_results_table(pscan_full_bg(my_gene_ids, full_bg))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.
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