## ----, echo = FALSE, results = 'hide'------------------------------------
library(knitr)
opts_chunk$set(error = FALSE)

## ----style, echo = FALSE, results = 'asis'-------------------------------
##BiocStyle::markdown()

## ----, message = FALSE---------------------------------------------------
library(SGSeq)

## ------------------------------------------------------------------------
si

## ------------------------------------------------------------------------
path <- system.file("extdata", package = "SGSeq")
si$file_bam <- file.path(path, "bams", si$file_bam)

## ----, message = FALSE---------------------------------------------------
library(TxDb.Hsapiens.UCSC.hg19.knownGene)
txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene
txdb <- keepSeqlevels(txdb, "chr16")
seqlevelsStyle(txdb) <- "NCBI"

## ------------------------------------------------------------------------
txf_ucsc <- convertToTxFeatures(txdb)
txf_ucsc <- txf_ucsc[txf_ucsc %over% gr]
txf_ucsc

## ------------------------------------------------------------------------
sgf_ucsc <- convertToSGFeatures(txf_ucsc)
sgf_ucsc

## ------------------------------------------------------------------------
sgfc_ucsc <- analyzeFeatures(si, features = txf_ucsc)
sgfc_ucsc

## ------------------------------------------------------------------------
head(counts(sgfc_ucsc))
head(FPKM(sgfc_ucsc))

## ----figure-1, fig.width=4.5, fig.height=4.5-----------------------------
df <- plotFeatures(sgfc_ucsc, geneID = 1)
df

## ------------------------------------------------------------------------
sgfc_pred <- analyzeFeatures(si, which = gr)

## ------------------------------------------------------------------------
sgfc_pred <- annotate(sgfc_pred, txf_ucsc)

## ----figure-2, fig.width=4.5, fig.height=4.5-----------------------------
df <- plotFeatures(sgfc_pred, geneID = 1, color_novel = "red")
df

## ------------------------------------------------------------------------
sgvc_pred <- analyzeVariants(sgfc_pred)
sgvc_pred

## ------------------------------------------------------------------------
variantFreq(sgvc_pred)

## ------------------------------------------------------------------------
mcols(sgvc_pred)

## ----figure-3, fig.width=1.5, fig.height=4.5-----------------------------
plotVariants(sgvc_pred, eventID = 1, color_novel = "red")

## ----, eval = FALSE------------------------------------------------------
## plotFeatures(sgfc_pred, geneID = 1)
## plotFeatures(sgfc_pred, geneName = "79791")
## plotFeatures(sgfc_pred, which = gr)

## ----, eval = FALSE------------------------------------------------------
## plotFeatures(sgfc_pred, geneID = 1, include = "junctions")
## plotFeatures(sgfc_pred, geneID = 1, include = "exons")
## plotFeatures(sgfc_pred, geneID = 1, include = "both")

## ----, eval = FALSE------------------------------------------------------
## plotFeatures(sgfc_pred, geneID = 1, toscale = "gene")
## plotFeatures(sgfc_pred, geneID = 1, toscale = "exon")
## plotFeatures(sgfc_pred, geneID = 1, toscale = "none")

## ----, figure-4, fig.width=4.5, fig.height=4.5---------------------------
par(mfrow = c(5, 1), mar = c(1, 3, 1, 1))
plotSpliceGraph(rowRanges(sgfc_pred), geneID = 1, toscale = "none", color_novel = "red")
for (j in 1:4) {
  plotCoverage(sgfc_pred[, j], geneID = 1, toscale = "none")
}

## ------------------------------------------------------------------------
txf <- predictTxFeatures(si, gr)
sgf <- convertToSGFeatures(txf)
sgf <- annotate(sgf, txf_ucsc)
sgfc <- getSGFeatureCounts(si, sgf)
sgv <- findSGVariants(sgf)
sgvc <- getSGVariantCounts(sgv, sgfc)

## ------------------------------------------------------------------------
sessionInfo()