## ----setup, include = FALSE--------------------------------------------------- options(tinytex.verbose = TRUE) knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) ## ----echo=FALSE, out.width=700------------------------------------------------ knitr::include_graphics('../README_files/VaSP.png') ## ----eval=FALSE--------------------------------------------------------------- # if (!requireNamespace("BiocManager", quietly = TRUE)) # install.packages("BiocManager") # BiocManager::install("VaSP") # vignette('VaSP') ## ----eval=FALSE--------------------------------------------------------------- # BiocManager::install("yuhuihui2011/VaSP", build_vignettes=TRUE) # vignette('VaSP') ## ----eval=FALSE--------------------------------------------------------------- # library(VaSP) # ?ballgown # path<-system.file('extdata', package='VaSP') # rice.bg<-ballgown(samples = list.dirs(path = path,recursive = F) ) ## ----------------------------------------------------------------------------- library(VaSP) data(rice.bg) ?rice.bg rice.bg score<-spliceGene(rice.bg, gene="MSTRG.183", junc.type = "score") tail(round(score,2),2) ## ----------------------------------------------------------------------------- gss <- BMfinder(score, cores = 1) gss ## ----------------------------------------------------------------------------- gss_intron<-structure(rice.bg)$intron (gss_intron<-gss_intron[gss_intron$id%in%rownames(gss)]) range(gss_intron) ## ----splicePlot, fig.width=7, fig.height=4, eval=FALSE------------------------ # splicePlot(rice.bg,gene='MSTRG.183',samples = sampleNames(rice.bg)[c(1,3,5)], # start = 1179000, end = 1179300) ## ----plotBedgraph, fig.height=3, fig.width=7---------------------------------- path <- system.file("extdata", package = "VaSP") bam_files <- list.files(path, "*.bam$") bam_files depth <- getDepth(file.path(path, bam_files[1]), "Chr1", start = 1171800, end = 1179400) head(depth) # library(Sushi) # par(mar=c(3,5,1,1)) # plotBedgraph(depth, "Chr1", chromstart = 1171800, chromend = 1179400,yaxt = "s") # mtext("Depth", side = 2, line = 2.5, cex = 1.2, font = 2) # labelgenome("Chr1", 1171800, 1179400, side = 1, scipen = 20, n = 5,scale = "Kb") ## ----plotGenes, fig.height=4,fig.width=7-------------------------------------- unique(geneIDs(rice.bg)) gene_id <- c("MSTRG.181", "MSTRG.182", "MSTRG.183") geneinfo <- getGeneinfo(genes = gene_id, rice.bg) trans <- table(geneinfo$name) # show how many exons each transcript has trans # chrom = geneinfo$chrom[1] # chromstart = min(geneinfo$start) - 1500 # chromend = max(geneinfo$stop) + 1000 # color = rep(SushiColors(2)(length(trans)), trans) # par(mar=c(3,1,1,1)) # p<-plotGenes(geneinfo, chrom, chromstart, chromend, col = color, bheight = 0.2, # bentline = FALSE, plotgenetype = "arrow", labeloffset = 0.5) # labelgenome(chrom, chromstart , chromend, side = 1, n = 5, scale = "Kb") ## ----------------------------------------------------------------------------- rice.bg head(geneIDs(rice.bg)) score <- spliceGene(rice.bg, "MSTRG.183", junc.type = "score") count <- spliceGene(rice.bg, "MSTRG.183", junc.type = "count") ## compare tail(score) tail(count) ## get intron structrue intron <- structure(rice.bg)$intron intron[intron$id %in% rownames(score)] ## ----------------------------------------------------------------------------- rice.bg splice <- spliceGenome(rice.bg, gene.select = NA, intron.select = NA) names(splice) head(splice$score) splice$intron ## ----------------------------------------------------------------------------- score <- spliceGene(rice.bg, "MSTRG.183", junc.type = "score") score <- round(score, 2) as <- BMfinder(score, cores = 1) # 4 bimodal distrubition features found ## compare as score[rownames(score) %in% rownames(as), ] ## ----fig.width=7, fig.height=4, eval=FALSE------------------------------------ # samples <- paste("Sample", c("027", "102", "237"), sep = "_") # bam.dir <- system.file("extdata", package = "VaSP") # # ## plot the whole gene region without junction lables # splicePlot(rice.bg, samples, bam.dir, gene = "MSTRG.183", junc.text = FALSE, # bheight = 0.2) # # ## plot the alternative splicing region with junction splicing scores # splicePlot(rice.bg, samples, bam.dir, gene = "MSTRG.183", start = 1179000) ## ----fig.width=7, fig.height=4, eval=FALSE------------------------------------ # splicePlot(rice.bg, samples, bam.dir, gene = "MSTRG.183", junc.type = 'count', # start = 1179000) ## ----------------------------------------------------------------------------- sessionInfo()