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