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

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

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

## ---- eval = FALSE-------------------------------------------------------
#  tx_1 <- GRangesList(tx_1 = GRanges("1", IRanges(c(101, 301, 501), c(200, 400, 600)), "+"))
#  tx_2 <- GRangesList(tx_2 = GRanges("1", IRanges(c(101, 351, 701), c(200, 400, 800)), "+"))
#  txf_1 <- convertToTxFeatures(tx_1)
#  txf_2 <- convertToTxFeatures(tx_2)
#  txf <- convertToTxFeatures(c(tx_1, tx_2))
#  sgf <- convertToSGFeatures(txf)
#  par(mfrow = c(1, 1))
#  plotSpliceGraph(sgf)

## ---- eval = FALSE-------------------------------------------------------
#  sgfc_IBM <- analyzeFeatures(si_IBM, alpha = 5, psi = 0.2, beta = 0.2, gamma = 0.2)
#  sgvc_IBM <- analyzeVariants(sgfc_IBM, min_denominator = 20)
#  exclude <- eventID(sgvc_IBM)[!closed5p(sgvc_IBM) | !closed3p(sgvc_IBM)]
#  sgvc_IBM <- sgvc_IBM[!eventID(sgvc_IBM) %in% exclude, ]

## ---- eval = FALSE-------------------------------------------------------
#  data(sgfc_IBM, sgvc_IBM, package = "SGSeqBioC2015")

## ---- eval = FALSE-------------------------------------------------------
#  txdb <- restoreSeqlevels(txdb)
#  seqlevelsStyle(txdb) <- "NCBI"
#  txdb <- keepSeqlevels(txdb, c(1:22, "X", "Y"))
#  txf_ucsc <- convertToTxFeatures(txdb)
#  sgfc_IBM <- annotate(sgfc_IBM, txf_ucsc)
#  sgvc_IBM <- annotate(sgvc_IBM, txf_ucsc)

## ---- eval = FALSE-------------------------------------------------------
#  plotFeatures(sgfc_IBM, geneName = "22920", color_novel = "red")
#  plotFeatures(sgfc_IBM, geneName = "22920", color_novel = "red", include = "exons")

## ---- eval = FALSE-------------------------------------------------------
#  mcols(sgvc_IBM)[any(geneName(sgvc_IBM) == "22920"), ]
#  plotVariants(sgvc_IBM, eventID = 552, color_novel = "red")

## ---- eval = FALSE-------------------------------------------------------
#  selected <- which(elementLengths(txName(sgvc_IBM)) == 0 & any(variantType(sgvc_IBM) == "SE:I"))
#  variants <- rowRanges(sgvc_IBM)[selected]
#  features <- unlist(variants)
#  exons <- features[type(features) == "E"]
#  exons_FPKM <- FPKM(sgfc_IBM)[match(featureID(exons), featureID(sgfc_IBM)), ]
#  exons_FPKM_max <- apply(exons_FPKM, 1, max)
#  geneIDs <- geneID(exons)[order(exons_FPKM_max, decreasing = TRUE)]
#  plotFeatures(sgfc_IBM, geneID = geneIDs[1], color_novel = "red")

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