--- title: "PscanR Concepts and Input Preparation" author: - Diana Betelli - Federico Zambelli date: "`r Sys.Date()`" output: BiocStyle::html_document: toc: true toc_depth: 2 vignette: > %\VignetteIndexEntry{2. PscanR Concepts and Input Preparation} %\VignetteEncoding{UTF-8} %\VignetteEngine{knitr::rmarkdown} editor_options: markdown: wrap: 80 --- ```{r setup, include = FALSE} library(PscanR) knitr::opts_chunk$set( collapse = TRUE, comment = "#>", warning = FALSE, message = FALSE ) ``` # 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 is an R implementation of the algorithm developed by [Zambelli et al. (2009)](https://doi.org/10.1093/nar/gkp464) for analysing sequence motif enrichment. It identifies transcription factors (TFs) that are potential regulators of a set of co-regulated or co-expressed genes. By scanning input sequences (e.g., promoters) against a database of known TF binding motifs, PscanR tests for motif enrichment and reports candidate binding sites. ## Relationship to other motif tools [TFBSTools](https://bioconductor.org/packages/TFBSTools/) supplies the motif classes that PscanR extends; Biostrings supplies sequence representations and scoring primitives. [motifmatchr](https://bioconductor.org/packages/motifmatchr/) uses MOODS to find matches across many motifs and sequences. [PWMEnrich](https://bioconductor.org/packages/PWMEnrich/) provides another Biostrings-based toolkit for known-motif scanning and enrichment. These tools address related tasks, but their statistical approach is different. PscanR implements the specific best-hit model of [Zambelli et al. (2009)](https://doi.org/10.1093/nar/gkp464): each sequence contributes its best normalized score per motif; the foreground mean is tested against promoter-background moments for a matching reference and window. Its `PSMatrix` and `PSMatrixList` classes extend TFBSTools objects to retain those moments and the selected positions, strands, and oligos. Gene-to-promoter selection and cached immutable backgrounds connect this model to workflows using annotation packages and `DNAStringSet` inputs. The full-background path also supports identifier-based queries without rescanning sequences. ## Workflow overview **Input:** DNA sequences supplied as a `DNAStringSet`, and transcription factor binding matrices with matching background statistics supplied as a `PSMatrixList`. **Analysis:** Calculation of the best motif match in each promoter, followed by comparison of foreground and background scores. **Output:** Motif-level z-scores, p-values, and adjusted p-values, together with the best match in each promoter for each motif. # The background model Statistical assessment requires a background model, represented by background statistics stored in a `PSMatrixList` object. PscanR treats the input sequences as a sample from a defined reference population, usually the annotated promoters of the species under investigation. This reference population supplies the background model for statistical analyses. Because the background is the distribution of scores over the whole promoter universe, the foreground and the background must describe the *same* kind of sequence. A promoter window of 200 bp upstream and 50 bp downstream of the transcription start site (TSS) does not yield scores comparable with one of 950 bp upstream and 50 bp downstream: a longer window offers more candidate matches and can therefore increase the best score without biological enrichment. Extract foreground promoters using the background's window and compatible reference sequences and annotations. Each precomputed background records its window in its name, as in `950u_50d`. The next section demonstrates the use of a precomputed background. The vignette [Building a background for PscanR](PscanR_Custom_Background.html) describes the construction and use of custom backgrounds. ## Obtaining a precomputed background PscanR provides access to precomputed backgrounds for five model organisms: *Homo sapiens*, *Mus musculus*, *Drosophila melanogaster*, *Arabidopsis thaliana*, and *Saccharomyces cerevisiae*. `ps_available_bg()` lists the available backgrounds, and `ps_retrieve_bg()` retrieves the selected background as a `PSMatrixList` for use with `pscan()`: ```{r list-backgrounds} ps_available_bg() ``` ```{r download-background} J2020_PSBG <- ps_retrieve_bg( "Jaspar2020", "hs", c(-200, 50), "hg38", source = "zenodo" ) J2020_PSBG J2020_PSBG[[4]] # one matrix, with its background statistics ``` The four arguments above identify the JASPAR release, species, promoter window relative to the TSS, and genome assembly. The assembly must be specified for mouse and human, for which multiple assemblies are available. The supported species and assemblies are: - *Homo sapiens*, assemblies hg38 and hs1 - *Mus musculus*, assemblies mm10 and mm39 - *Arabidopsis thaliana*, TAIR9 - *Drosophila melanogaster*, dm6 - *Saccharomyces cerevisiae*, sacCer3 The available promoter windows span 200 bp upstream and 50 bp downstream, 450 upstream and 50 downstream, 500 upstream and 0 downstream, 950 upstream and 50 downstream, and 1000 upstream and 0 downstream, for the JASPAR 2020, 2022, and 2024 releases. Only JASPAR CORE motifs are used. Each precomputed background is a separate ExperimentHub resource of the `PscanRBackgrounds` data package, and the files are hosted on Zenodo in one record per JASPAR release. By default `ps_retrieve_bg()` retrieves a background through ExperimentHub and, if the Hub cannot be reached, falls back to Zenodo with a warning; `source = "zenodo"`, used above, downloads the file from Zenodo directly. If a single file cannot be retrieved, the immutable record , which archives all backgrounds in a single ZIP, is used instead. Downloads are cached, and every file is verified against the SHA-256 checksum recorded in the catalog. For offline use, the package also includes JASPAR 2020 CORE matrices and a historical human background table for 200 bp upstream and 50 bp downstream of the TSS, computed on an earlier snapshot of the hg38 annotation. `ps_retrieve_bg_from_file()` combines them into a `PSMatrixList`; the man page examples use it: ```{r offline-background} J2020 <- readRDS(system.file("extdata", "J2020.rds", package = "PscanR")) bg_path <- system.file( "extdata", "J2020_hg38_200u_50d_UCSC.psbg.txt", package = "PscanR" ) J2020_offline <- ps_retrieve_bg_from_file(bg_path, J2020) length(J2020_offline) ``` If no available background matches the required species, window, or motif collection, construct a custom background as described in [Building a background for PscanR](PscanR_Custom_Background.html). # Preparing promoter sequences for `pscan()` The two principal inputs to `pscan()` are a `DNAStringSet` of promoter sequences and a background `PSMatrixList`. The following sections demonstrate promoter selection from transcript identifiers and from gene identifiers. For humans, mice, and Drosophila, the available background models are based on RefSeq identifiers; for *Arabidopsis*, on TAIR; and for yeast, on SGD. Other types of gene identifiers can also be mapped to the annotation used by an existing background. Custom backgrounds can use other annotations, as described in [Building a background for PscanR](PscanR_Custom_Background.html). ## Scanning from a list of transcripts PscanR operates on promoter sequences. When transcript identifiers map uniquely within the chosen assembly and annotation, their corresponding promoters can be selected directly without choosing among transcripts of the same gene. The package includes the historical human NRF1 example list of RefSeq transcripts from the original Pscan web interface. See [Zambelli, Pesole and Pavesi (2009)](https://doi.org/10.1093/nar/gkp464) for the software and interface. ```{r chunk-2} file_path <- system.file("extdata", "nrf1100.txt", package = "PscanR") target <- read.csv(file_path, header = FALSE) head(target) ``` The bundled human promoter sequences span 200 bp upstream and 50 bp downstream of the TSS, matching the background loaded above: ```{r load-promoters} library(Biostrings) prom_seq <- readRDS( system.file("extdata", "prom_seq.rds", package = "PscanR") ) length(prom_seq) unique(width(prom_seq)) # 200 upstream + 50 downstream head(names(prom_seq), 3) ``` The transcript list and promoter names use identical RefSeq identifiers, including version suffixes. The foreground can therefore be selected by exact name matching: ```{r match-by-name} foreground <- prom_seq[names(prom_seq) %in% target[[1]]] length(foreground) ``` Every transcript in the list has a promoter in the bundled set, so all 90 are retained. The same operation can select a foreground from a genome-wide promoter set. The selected sequences are then scanned against the background: ```{r run-pscan} res <- pscan(foreground, J2020_PSBG, BPPARAM = BiocParallel::SerialParam()) head(ps_results_table(res)) ``` `ps_results_table()` returns motifs ordered by increasing p-value, with decreasing z-score used to resolve ties. NRF1 is the highest-ranked motif, consistent with the origin of the curated target set. ## Selecting promoters from a gene list A differential expression analysis often produces a gene list. A gene may have several transcription start sites and therefore several distinct promoters, so gene identifiers alone do not uniquely specify a promoter set. This workflow selects one promoter per gene to avoid giving greater weight to genes with more annotated transcripts. `ps_select_promoters()` applies an explicit selection rule and records the selected transcript and the reason for its selection. This provides reproducibility without establishing which transcript is active in the experimental condition. The function requires the genes of interest, candidate promoter sequences, and a gene-to-transcript annotation table. The following subset of genes from the NRF1 target list illustrates the input format: ```{r gene-list} genes <- c( "ALAS1", "ARPC1A", "CAPNS1", "CHMP2A", "EIF2S1", "GTF2F1", "LONP1", "POLR2A" ) ``` `org.Hs.eg.db` supplies the gene-to-transcript mapping. These eight genes have more than eight candidate transcripts: ```{r build-annotation} library(org.Hs.eg.db) annotation <- AnnotationDbi::select( org.Hs.eg.db, keys = genes, keytype = "SYMBOL", columns = "REFSEQ" ) annotation <- annotation[!is.na(annotation$REFSEQ), ] nrow(annotation) table(annotation$SYMBOL) ``` The annotation table contains unversioned RefSeq accessions, whereas promoter names such as `NM_006559.3` include a version suffix. The function recognises the RefSeq convention and removes the suffix for identifier matching. It preserves the suffix in TAIR identifiers such as `AT1G01110.2`, where it distinguishes transcript variants that may have different start sites. Thus, RefSeq matching uses the accession after version normalisation, whereas TAIR matching uses the complete transcript identifier. Matching RefSeq versions does not imply that their sequences or annotated start sites are identical; the promoter sequences must still correspond to the chosen reference. ```{r select-promoters} selection <- ps_select_promoters( genes = genes, promoter_sequences = prom_seq, annotation = annotation, gene_col = "SYMBOL", transcript_col = "REFSEQ", mode = "representative", return = "both" ) ``` The returned mapping documents each selection. Protein accessions such as `NP_` are excluded because they do not identify transcripts: ```{r selection-mapping} selection$mapping[, c( "gene", "transcript_id", "promoter_id", "scheme", "suffix_role", "n_candidates", "decided_by" )] ``` `ps_selection_summary()` summarises the mapping, including the number of genes with multiple candidates and the number requiring a tie-breaking rule after transcript ranking: ```{r selection-summary} ps_selection_summary(selection$mapping) ``` ```{r selection-counts, include=FALSE} unmapped <- attr(selection$mapping, "ps_unmapped_genes") n_selected <- length(selection$sequences) ``` The number of unmapped genes is `r length(unmapped)` out of `r length(genes)`; the unmapped identifiers are `r paste(unmapped, collapse = ", ")`. These genes have no matching transcript promoter in the bundled set. The `r n_selected` selected promoters are returned as a `DNAStringSet`, named by gene and transcript so that a result can be traced back to both: ```{r selected-sequences} selection$sequences ``` `mode = "representative"` ranks transcripts by the identifier scheme's rule. For RefSeq, it prefers `NM_` over `NR_`, followed by predicted `XM_` and `XR_` accessions. For human hg38, MANE Select and RefSeq Select provide curated transcript preferences. These can be used with `mode = "select"`: ```{r mane-select, eval=FALSE} select_tx <- ps_load_select_transcripts("hg38") selection <- ps_select_promoters( genes = genes, promoter_sequences = prom_seq, annotation = annotation, select_transcripts = select_tx, mode = "select", return = "both" ) ``` `ps_load_select_transcripts()` retrieves the two tracks through the UCSC API. This network-dependent example is not evaluated during vignette compilation. After retrieval, the table is cached under `tools::R_user_dir("PscanR", "cache")` for subsequent calls. The metadata is specific to human hg38, so for other organisms either `mode = "representative"` applies, or a species-specific table is supplied through `select_transcripts`. With `return = "both"`, the selected `DNAStringSet` is available as `selection$sequences` and can be passed to `pscan()`. It contains one promoter per mapped gene. The sequences and background must use the same window, here 200 bp upstream and 50 bp downstream. A `950u_50d` background would require promoters extracted over the corresponding interval. # Configuring parallel execution `pscan()` parallelises over motifs through **BiocParallel**, and the `BPPARAM` argument selects the execution strategy: - **`MulticoreParam(n)` on systems supporting forked execution.** Workers inherit the parent process's memory, which can reduce data-transfer overhead on Linux and macOS. Suitability also depends on the execution environment; forked execution is unavailable on Windows. - **`SnowParam(n)` for portable parallel execution.** Independent R processes communicate through sockets on Windows and Unix-like systems. Worker startup and data transfer introduce additional overhead. - **`SerialParam()` for debugging and small examples.** Computation occurs within the current R process, simplifying error inspection. The default is `BiocParallel::SerialParam()`, which uses one R process. A registered backend does not change this default; pass `BPPARAM` explicitly to request multiple workers. Select the worker count according to the workload and available resources. For short scans, worker startup and data-transfer costs may exceed the benefit of parallel execution. # 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() ```