## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
    collapse = TRUE,
    comment = "#>"
)

## ----library, results='hide', message=FALSE-----------------------------------
library(TSSr)

## ----citation, eval=TRUE------------------------------------------------------
citation("TSSr")

## ----bundled-input-files------------------------------------------------------
exampleInput <- system.file(
    "extdata",
    "example-tss-table.tsv",
    package = "TSSr",
    mustWork = TRUE
)
exampleAnnotation <- system.file(
    "extdata",
    "example-annotation.gff3",
    package = "TSSr",
    mustWork = TRUE
)

## ----construct-object---------------------------------------------------------
myTSSr <- TSSr(
    genomeName = "BSgenome.Scerevisiae.UCSC.sacCer3",
    inputFiles = exampleInput,
    inputFilesType = "TSStable",
    sampleLabels = c("SL01", "SL02", "SL03", "SL04"),
    sampleLabelsMerged = c("control", "treat"),
    mergeIndex = c(1, 1, 2, 2),
    refSource = exampleAnnotation
)

myTSSr

## ----import-tss, message=FALSE------------------------------------------------
myTSSr <- getTSS(myTSSr)
dim(TSSmatrix(myTSSr, data = "raw"))

## ----TSSprocessing, tidy=FALSE------------------------------------------------
# Merge replicates
myTSSr <- mergeSamples(myTSSr)
# Normalization
myTSSr <- normalizeTSS(myTSSr)
# TSS filtering
rowsBeforeFiltering <- nrow(TSSmatrix(myTSSr, data = "processed"))
myTSSr <- filterTSS(myTSSr, method = "TPM", tpmLow = 2)
c(
    before = rowsBeforeFiltering,
    after = nrow(TSSmatrix(myTSSr, data = "processed"))
)

## ----TSSprocessing-show-------------------------------------------------------
# Access the processed TSS matrix without exposing internal storage
head(TSSmatrix(myTSSr, data = "processed"))

## ----TSSclustering, tidy=FALSE------------------------------------------------
# TSS clustering
myTSSr <- clusterTSS(myTSSr,
    method = "peakclu", peakDistance = 100, extensionDistance = 30,
    localThreshold = 0.02, clusterThreshold = 1,
    useMultiCore = FALSE, numCores = NULL
)

# Aggregating consensus clusters
myTSSr <- consensusCluster(myTSSr, dis = 50, useMultiCore = FALSE)

## ----TSSclustering-show-------------------------------------------------------
# Tag clusters per sample
head(tagClusters(myTSSr, sample = "control"))

# Consensus clusters across samples
head(consensusClusters(myTSSr, sample = "control"))

## ----shapeCluster, tidy=FALSE-------------------------------------------------
# Calculating core promoter shape score
myTSSr <- shapeCluster(myTSSr,
    clusters = "consensusClusters", method = "PSS",
    useMultiCore = FALSE, numCores = NULL
)

## ----shapeCluster-show--------------------------------------------------------
# Shape scores per cluster
head(clusterShape(myTSSr, sample = "control"))

## ----annotateCluster, tidy=FALSE----------------------------------------------
# Assign clusters to the annotated features
myTSSr <- annotateCluster(myTSSr, clusters = "consensusClusters",
    filterCluster = TRUE, filterClusterThreshold = 0.02,
    annotationType = "genes", upstream = 1000,
    upstreamOverlap = 500, downstream = 0)

## ----annotateCluster-show-----------------------------------------------------
# Clusters assigned to genes
head(assignedClusters(myTSSr, sample = "control"))

## ----callEnhancer, tidy=FALSE-------------------------------------------------
myTSSr <- callEnhancer(myTSSr, flanking = 400, dis2gene = 2000)

## ----callEnhancer-show--------------------------------------------------------
# Putative enhancers per sample
head(enhancers(myTSSr, sample = "control"))

## ----deGene, tidy=FALSE-------------------------------------------------------
# Gene-level differential expression using DESeq2
myTSSr <- deGene(myTSSr, comparePairs = list(c("control", "treat")),
    pval = 0.01, useMultiCore = FALSE, numCores = NULL,
    fitType = "mean")

## ----deGene-show--------------------------------------------------------------
# All differential expression results
head(DEtables(
    myTSSr,
    comparison = "control_VS_treat",
    result = "all"
))

## ----shiftPromoter, tidy=FALSE------------------------------------------------
# Calculate core promoter shifts
myTSSr <- shiftPromoter(myTSSr, comparePairs = list(c("control", "treat")),
    pval = 0.01)

## ----shiftPromoter-show-------------------------------------------------------
# Promoter shift results
head(PromoterShift(myTSSr, comparison = "control_VS_treat"))

## ----export-------------------------------------------------------------------
oldwd <- setwd(tempdir())

# Export processed TSS values to ALL.samples.TSS.processed.txt
exportTSStable(myTSSr, data = "processed")

# Export cluster tables
exportClustersTable(myTSSr, data = "assigned")

# Export to BED format for genome browsers
exportClustersToBed(myTSSr, data = "consensusClusters")

# Export TSS to bedGraph
exportTSStoBedgraph(myTSSr, data = "processed")

setwd(oldwd)

## ----precomputed-example------------------------------------------------------
data("exampleTSSr")
exampleTSSr

## ----complete-pipeline, eval=FALSE--------------------------------------------
# myTSSr <- TSSr(
#     genomeName = "BSgenome.Scerevisiae.UCSC.sacCer3",
#     inputFiles = exampleInput,
#     inputFilesType = "TSStable",
#     sampleLabels = c("SL01", "SL02", "SL03", "SL04"),
#     sampleLabelsMerged = c("control", "treat"),
#     mergeIndex = c(1, 1, 2, 2),
#     refSource = exampleAnnotation
# ) |>
#     getTSS() |>
#     mergeSamples() |>
#     normalizeTSS() |>
#     filterTSS(method = "TPM", tpmLow = 2) |>
#     clusterTSS(
#         method = "peakclu",
#         peakDistance = 100,
#         extensionDistance = 30,
#         localThreshold = 0.02,
#         clusterThreshold = 1,
#         useMultiCore = FALSE,
#         numCores = NULL
#     ) |>
#     consensusCluster(
#         dis = 50,
#         useMultiCore = FALSE
#     ) |>
#     shapeCluster(
#         clusters = "consensusClusters",
#         method = "PSS",
#         useMultiCore = FALSE,
#         numCores = NULL
#     ) |>
#     annotateCluster(
#         clusters = "consensusClusters",
#         filterCluster = TRUE,
#         filterClusterThreshold = 0.02,
#         annotationType = "genes",
#         upstream = 1000,
#         upstreamOverlap = 500,
#         downstream = 0
#     ) |>
#     callEnhancer(
#         flanking = 400,
#         dis2gene = 2000
#     ) |>
#     deGene(
#         comparePairs = list(c("control", "treat")),
#         pval = 0.01,
#         useMultiCore = FALSE,
#         numCores = NULL,
#         fitType = "mean"
#     ) |>
#     shiftPromoter(
#         comparePairs = list(c("control", "treat")),
#         pval = 0.01
#     )

## ----sessionInfo, tidy=FALSE--------------------------------------------------
sessionInfo()