## ----setup, include = FALSE---------------------------------------------------
library(PscanR)

knitr::opts_chunk$set(
    collapse = TRUE,
    comment = "#>",
    warning = FALSE,
    message = FALSE
)

## ----installation, eval=FALSE-------------------------------------------------
# BiocManager::install("PscanR")

## ----list-backgrounds---------------------------------------------------------
ps_available_bg()


## ----download-background------------------------------------------------------
J2020_PSBG <- ps_retrieve_bg(
    "Jaspar2020", "hs", c(-200, 50), "hg38"
)
J2020_PSBG
J2020_PSBG[[4]] # one matrix, with its background statistics

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

## ----chunk-2------------------------------------------------------------------
file_path <- system.file("extdata", "nrf1100.txt", package = "PscanR")
target <- read.csv(file_path, header = FALSE)
head(target)

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

## ----match-by-name------------------------------------------------------------
foreground <- prom_seq[names(prom_seq) %in% target[[1]]]
length(foreground)

## ----run-pscan----------------------------------------------------------------
res <- pscan(foreground, J2020_PSBG, BPPARAM = BiocParallel::SerialParam())
head(ps_results_table(res))

## ----gene-list----------------------------------------------------------------
genes <- c(
    "ALAS1", "ARPC1A", "CAPNS1", "CHMP2A",
    "EIF2S1", "GTF2F1", "LONP1", "POLR2A"
)

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

## ----select-promoters---------------------------------------------------------
selection <- ps_select_promoters(
    genes = genes,
    promoter_sequences = prom_seq,
    annotation = annotation,
    gene_col = "SYMBOL", transcript_col = "REFSEQ",
    mode = "representative",
    return = "both"
)

## ----selection-mapping--------------------------------------------------------
selection$mapping[, c(
    "gene", "transcript_id", "promoter_id",
    "scheme", "suffix_role", "n_candidates", "decided_by"
)]

## ----selection-summary--------------------------------------------------------
ps_selection_summary(selection$mapping)

## ----selection-counts, include=FALSE------------------------------------------
unmapped <- attr(selection$mapping, "ps_unmapped_genes")
n_selected <- length(selection$sequences)

## ----selected-sequences-------------------------------------------------------
selection$sequences

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

## ----session-info-------------------------------------------------------------
sessionInfo()