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

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


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

## ----read-gene-sets-----------------------------------------------------------
membership <- read.delim(
    system.file(
        "extdata", "vignettes", "mouse_hd_retina",
        "paper_subset_membership.tsv",
        package = "PscanR"
    ),
    stringsAsFactors = FALSE
)
str(membership)
gene_sets <- split(membership$ensembl_gene_id, membership$subset)
lengths(gene_sets)

## ----subset-key, echo=FALSE---------------------------------------------------
subset_labels <- c(
    a = "retina-only down", b = "shared down", c = "striatum-only down",
    d = "retina-only up", e = "shared up", f = "striatum-only up"
)
knitr::kable(
    data.frame(
        set = names(subset_labels),
        definition = unname(subset_labels),
        genes = lengths(gene_sets)[names(subset_labels)],
        row.names = NULL
    ),
    caption = "The six gene sets, by tissue and direction of change."
)

## ----discover-background------------------------------------------------------
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."
)

## ----retrieve-background------------------------------------------------------
background <- ps_retrieve_bg(
    jaspar_release = "Jaspar2020",
    org = "mm",
    prom_reg = c(-950, 50),
    assembly = "mm10",
    version = 2
)
length(background)

## ----promoter-scheme, message=TRUE--------------------------------------------
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"
)

## ----promoter-scheme-table, echo=FALSE----------------------------------------
knitr::kable(
    selection[, c(
        "gene", "promoter_id", "scheme", "suffix_role",
        "transcript_class", "selection_source", "n_candidates", "decided_by"
    )],
    row.names = FALSE,
    caption = "Selected promoters and the rules used for their selection."
)

## ----select-real, eval=FALSE--------------------------------------------------
# 
# 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]

## ----load-precomputed---------------------------------------------------------
precomputed <- readRDS(system.file(
    "extdata", "vignettes", "mouse_hd_retina", "uc2_vignette_data.rds",
    package = "PscanR"
))
retina_mapping <- precomputed$retina_mapping

## ----selection-summary--------------------------------------------------------
ps_selection_summary(retina_mapping)$overall

## ----unmapped-----------------------------------------------------------------
head(attr(retina_mapping, "ps_unmapped_genes"), 4)

## ----selection-by-source------------------------------------------------------
ps_selection_summary(retina_mapping)$by_source

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

## ----scan-time, include=FALSE-------------------------------------------------
timings <- precomputed$timings
scan_seconds <- round(
    sum(timings$elapsed_seconds[timings$subset %in% c("d", "f")]), 0
)

## ----scan, eval=FALSE---------------------------------------------------------
# 
# 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
# )

## ----results-table, eval=FALSE------------------------------------------------
# enrichment <- ps_results_table(results)

## ----load-results-------------------------------------------------------------
# Saved enrichment table for retina-only upregulated promoters.
enrichment <- precomputed$enrichment_d
dim(enrichment)

# A motif subset of the PSMatrixList retains the per-promoter results.
results_subset <- precomputed$results_subset_d
results_subset

## ----subset-names-------------------------------------------------------------
TFBSTools::name(results_subset)

## ----results-table-shown------------------------------------------------------
head(enrichment[order(-enrichment$ZSCORE), ], 8)

## ----bg-table-----------------------------------------------------------------
ps_get_bg_table(results_subset)

## ----top-motifs, fig.width=8, fig.height=5, dpi=120---------------------------

table(ps_motif_class(results_subset))

ps_motif_barplot(results_subset, group = "class")

## ----barplot-composable-------------------------------------------------------
ps_motif_barplot(results_subset, group = "class") +
    ggplot2::labs(title = "Retina-only upregulated promoters")

## ----second-set-table---------------------------------------------------------
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))
)

## ----second-set-barplot, fig.width=8, fig.height=5, dpi=120-------------------
ps_motif_barplot(precomputed$results_subset_f, group = "class")

## ----irf-stat-data------------------------------------------------------------
irf_stat <- precomputed$irf_stat
str(irf_stat[c("subset", "motif", "ZSCORE", "significance")])

## ----irf-stat-heatmap, fig.width=10, fig.height=6.7, dpi=120------------------
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"
    )

## ----irf-stat-counts, echo=FALSE----------------------------------------------
counts <- do.call(rbind, lapply(split(irf_stat, irf_stat$subset), function(x) {
    data.frame(
        set = as.character(x$subset[[1]]),
        `nominal P < 0.05` = sum(x$P.VALUE < 0.05),
        `FDR < 0.05` = sum(x$FDR < 0.05),
        `mean Z` = round(mean(x$ZSCORE), 2),
        check.names = FALSE
    )
}))
knitr::kable(
    counts, row.names = FALSE,
    caption = paste(
        "The eighteen IRF and STAT matrices, counted in each set."
    )
)

## ----subset-e-size, include=FALSE---------------------------------------------
qc <- precomputed$subset_qc
n_subset_e <- qc$scan_eligible_sequences[qc$subset == "e"]

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

ps_zscore(irf3)
ps_pvalue(irf3)
c(background = ps_bg_avg(irf3), foreground = ps_fg_avg(irf3))
ps_hits_size(irf3)

## ----hits-table---------------------------------------------------------------
head(ps_hits_table(irf3, pos_shift = -950L), 5)

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

## ----z-table------------------------------------------------------------------
z <- ps_z_table(results_subset)
dim(z)
round(z[1:4, 1:5], 2)

## ----zscore-heatmap, fig.width=9, fig.height=7, dpi=110-----------------------
ps_zscore_heatmap(
    results_subset, FDR = 1,
    show_rownames = FALSE, treeheight_row = 0
)

## ----hitpos-map, fig.width=9, fig.height=7, dpi=110---------------------------
ps_hitpos_map(
    results_subset, FDR = 1, shift = -950,
    show_rownames = FALSE, treeheight_row = 0
)

## ----density-plot, fig.width=7.5, fig.height=4.6, dpi=120---------------------
ps_density_plot(
    irf3,
    shift = -950, st = "strict",
    window = c(-950, 50), bins = 20
)

## ----density-counts, include=FALSE--------------------------------------------
irf3_scores <- ps_hits_score(irf3)
loose_n <- sum(irf3_scores >= ps_bg_avg(irf3))
strict_n <- sum(irf3_scores >= ps_bg_avg(irf3) + ps_bg_std_dev(irf3))

## ----hit-score-plot, fig.width=9.5, fig.height=4.6, dpi=115-------------------
ps_hit_score_plot(
    results_subset[c("MA1418.1", "MA1513.1")],
    shift = -950
)

## ----pscan-filtered, eval=FALSE-----------------------------------------------
# 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)

## ----filtered-table, echo=FALSE-----------------------------------------------
ft <- precomputed$filtered_irf3$table
ft <- head(ft[order(-ft$ZSCORE), ], 12)
knitr::kable(
    data.frame(
        motif_id = ft$motif_id, NAME = ft$NAME,
        ZSCORE = round(ft$ZSCORE, 2),
        P.VALUE = format(ft$P.VALUE, scientific = TRUE, digits = 3)
    ),
    row.names = FALSE, align = c("l", "l", "r", "r"),
    caption = paste(
        "Enrichment within the", precomputed$filtered_irf3$n_retained,
        "promoters with a high-scoring IRF3 match."
    )
)

## ----filtered-contrast, echo=FALSE--------------------------------------------
zinc <- c("MA1513.1", "MA1650.1", "MA0162.4")
full <- precomputed$filtered_irf3$table
knitr::kable(
    data.frame(
        NAME = enrichment[zinc, "NAME"],
        `Z all promoters` = round(enrichment[zinc, "ZSCORE"], 2),
        `Z IRF3-positive` = round(
            full$ZSCORE[match(zinc, full$motif_id)], 2
        ),
        check.names = FALSE
    ),
    row.names = FALSE,
    caption = paste(
        "GC-rich motif scores before and after promoter selection",
        "by IRF3 score."
    )
)

## ----workflow, eval=FALSE-----------------------------------------------------
# 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)

## ----selection-counts, include=FALSE------------------------------------------
retina_selection <- ps_selection_summary(retina_mapping)$overall

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