inst/doc/SGSeq.R

In SGSeq: Splice event prediction and quantification from RNA-seq data

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