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