## ---- 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) ## ---- eval = FALSE------------------------------------------------------- # tx_1 <- GRangesList(tx_1 = GRanges("1", IRanges(c(101, 301, 501), c(200, 400, 600)), "+")) # txf_1 <- convertToTxFeatures(tx_1) # tx_2 <- GRangesList(tx_2 = GRanges("1", IRanges(c(101, 351, 701), c(200, 400, 800)), "+")) # 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 <- analyzeFeatures(si, alpha = 2, beta = 0.2, gamma = 0.2) # sgvc <- analyzeVariants(sgfc, min_denominator = 10) ## ---- eval = TRUE-------------------------------------------------------- data(sgfc, sgvc, package = "SGSeqBioC2016") txdb <- restoreSeqlevels(txdb) seqlevelsStyle(txdb) <- "NCBI" txdb <- keepSeqlevels(txdb, c(1:22, "X", "Y")) ## ---- eval = FALSE------------------------------------------------------- # sgf <- rowRanges(sgfc) # plotSpliceGraph(sgf, geneName = "22920") # plotSpliceGraph(sgf, geneName = "22920", color_novel = "red") # plotFeatures(sgfc, geneName = "22920", color_novel = "red") # plotFeatures(sgfc, geneName = "22920", color_novel = "red", include = "exons") ## ---- eval = FALSE------------------------------------------------------- # mcols(sgvc)[any(geneName(sgvc) == "22920"), ] # plotVariants(sgvc, eventID = 2872, color_novel = "red") ## ---- eval = FALSE------------------------------------------------------- # variant <- rowRanges(sgvc)[variantID(sgvc) == 6796] # variant <- keepSeqlevels(variant, c(1:22, "X", "Y")) # vep <- predictVariantEffects(variant, txdb, Hsapiens) # vep ## ------------------------------------------------------------------------ sessionInfo()