## ---- 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]
head(txf_ucsc)

## --------------------------------------------------------------------------
type(txf_ucsc)
head(txName(txf_ucsc))
head(geneName(txf_ucsc))

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

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

## ---- eval = FALSE---------------------------------------------------------
#  colData(sgfc_ucsc)
#  rowRanges(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)

## ---- message = FALSE------------------------------------------------------
sgfc_pred <- analyzeFeatures(si, which = gr)
head(rowRanges(sgfc_pred))

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

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

## ---- message = FALSE------------------------------------------------------
sgvc_pred <- analyzeVariants(sgfc_pred)
sgvc_pred

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

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

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

## ---- message = FALSE------------------------------------------------------
library(BSgenome.Hsapiens.UCSC.hg19)
seqlevelsStyle(Hsapiens) <- "NCBI"
vep <- predictVariantEffects(sgv_pred, txdb, Hsapiens)
vep

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

## ---- message = FALSE------------------------------------------------------
sgv <- rowRanges(sgvc_pred)
sgvc <- getSGVariantCounts(sgv, sample_info = si)
sgvc

## --------------------------------------------------------------------------
x <- counts(sgvc)
vid <- variantID(sgvc)
eid <- eventID(sgvc)

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

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