## ----eval=FALSE------------------------------------------------------------ # if (!"BiocManager" %in% rownames(installed.packages())) # install.packages("BiocManager") # BiocManager::install("APAlyzer") ## ----eval=FALSE------------------------------------------------------------ # R CMD INSTALL APAlyzer.tar.gz ## ----eval=FALSE------------------------------------------------------------ # BiocManager::install('RJWANGbioinfo/APAlyzer') ## -------------------------------------------------------------------------- library(APAlyzer) ## ----eval=TRUE------------------------------------------------------------- suppressMessages(library("TBX20BamSubset")) suppressMessages(library("Rsamtools")) flsall = getBamFileList() flsall ## ----eval=TRUE------------------------------------------------------------- extpath = system.file("extdata", "mm9_REF.RData", package="APAlyzer") load(extpath, verbose=TRUE) ## ----eval=TRUE------------------------------------------------------------- head(refUTRraw,2) ## ----eval=TRUE------------------------------------------------------------- head(dfIPA,2) ## ----eval=TRUE------------------------------------------------------------- head(dfLE,2) ## ----eval=TRUE------------------------------------------------------------- extpath = system.file("extdata", "hg19_REF.RData", package="APAlyzer") load(extpath, verbose=TRUE) ## ----eval=FALSE------------------------------------------------------------ # refUTRraw=refUTRraw[which(refUTRraw$Chrom=='chr19'),] # UTRdbraw=REF3UTR(refUTRraw) ## ---- echo = FALSE--------------------------------------------------------- extpath = system.file("extdata", "mm9_REF.RData", package="APAlyzer") load(extpath) refUTRraw=refUTRraw[which(refUTRraw$Chrom=='chr19'),] UTRdbraw=REF3UTR(refUTRraw) ## -------------------------------------------------------------------------- head(UTRdbraw,2) ## ---- echo = FALSE--------------------------------------------------------- options(warn=-1) suppressMessages(library("TBX20BamSubset")) suppressMessages(library("Rsamtools")) extpath = system.file("extdata", "mm9_REF.RData", package="APAlyzer") load(extpath) flsall <- getBamFileList() refUTRraw=refUTRraw[which(refUTRraw$Chrom=='chr19'),] UTRdbraw=REF3UTR(refUTRraw) ## ----eval=TRUE------------------------------------------------------------- DFUTRraw=PASEXP_3UTR(UTRdbraw, flsall, Strandtype="forward") ## ----eval=TRUE------------------------------------------------------------- head(DFUTRraw,2) ## ----eval=TRUE------------------------------------------------------------- #mouse(mm9): extpath = system.file("extdata", "mm9_REF.RData", package="APAlyzer") load(extpath) ## ----eval=TRUE------------------------------------------------------------- #human(hg19): extpath = system.file("extdata", "hg19_REF.RData", package="APAlyzer") load(extpath) ## ----eval=FALSE------------------------------------------------------------ # dfIPA=dfIPA[which(dfIPA$Chrom=='chr19'),] # dfLE=dfLE[which(dfLE$Chrom=='chr19'),] # IPA_OUTraw=PASEXP_IPA(dfIPA, dfLE, flsall, Strandtype="forward", nts=1) ## ---- echo = FALSE--------------------------------------------------------- extpath = system.file("extdata", "mm9_TBX20.APAout.RData", package="APAlyzer") load(extpath) ## -------------------------------------------------------------------------- head(IPA_OUTraw,2) ## ---- echo = FALSE--------------------------------------------------------- options(warn=-1) extpath = system.file("extdata", "mm9_TBX20.APAout.RData", package="APAlyzer") load(extpath) ## ----eval=TRUE------------------------------------------------------------- # Build the sample table with replicates sampleTable1 = data.frame(samplename = c(names(flsall)), condition = c(rep("NT",3),rep("KD",3))) sampleTable1 ## ----eval=TRUE------------------------------------------------------------- # Build the sample table without replicates sampleTable2 = data.frame(samplename = c("SRR316184","SRR316187"), condition = c("NT","KD")) sampleTable2 ## ----eval=TRUE------------------------------------------------------------- # Analysis 3'UTR APA between KD and NT group using non-repilicate design test_3UTRsing=APAdiff(sampleTable2,DFUTRraw, conKET='NT', trtKEY='KD', PAS='3UTR', CUTreads=0) ## ----eval=TRUE------------------------------------------------------------- head(test_3UTRsing,2) table(test_3UTRsing$APAreg) ## ----eval=TRUE------------------------------------------------------------- # Analysis 3'UTR APA between KD and NT group using multi-repilicate design test_3UTRmuti=APAdiff(sampleTable1, DFUTRraw, conKET='NT', trtKEY='KD', PAS='3UTR', CUTreads=0) head(test_3UTRmuti,2) table(test_3UTRmuti$APAreg) ## ----eval=TRUE------------------------------------------------------------- test_IPAsing=APAdiff(sampleTable2, IPA_OUTraw, conKET='NT', trtKEY='KD', PAS='IPA', CUTreads=0) head(test_IPAsing,2) ## ----eval=TRUE------------------------------------------------------------- test_IPAmuti=APAdiff(sampleTable1, IPA_OUTraw, conKET='NT', trtKEY='KD', PAS='IPA', CUTreads=0) head(test_IPAmuti,2) ## ---- echo = FALSE--------------------------------------------------------- options(warn=-1) extpath = system.file("extdata", "mm9_TBX20.APAdiff_OUT.RData", package="APAlyzer") load(extpath) ## ----eval=TRUE------------------------------------------------------------- test_3UTRmuti$APA="3'UTR" test_IPAmuti$APA="IPA" dfplot=rbind(test_3UTRmuti[,c('RED','APA')],test_IPAmuti[,c('RED','APA')]) ## ----eval=TRUE------------------------------------------------------------- library(ggplot2) ## ----fig1, fig.height = 4, fig.width = 4, fig.align = "center"------------- ###violin ggplot(dfplot, aes(x = APA, y = RED)) + geom_violin(trim = FALSE) + geom_boxplot(width = 0.2)+ theme_bw() + geom_hline(yintercept=0, linetype="dashed", color = "red") ## ----fig2, fig.height = 4, fig.width = 5, fig.align = "center"------------- ###CDF ggplot(dfplot, aes( x = RED, color = APA)) + stat_ecdf(geom = "step") + ylab("cumulative fraction")+ geom_vline(xintercept=0, linetype="dashed", color = "gray")+ theme_bw() + geom_hline(yintercept=0.5, linetype="dashed", color = "gray") ## ----eval=TRUE------------------------------------------------------------- suppressMessages(library("GenomicFeatures")) suppressMessages(library("org.Mm.eg.db")) extpath = system.file("extdata", "mm9.chr19.refGene.R.DB", package="APAlyzer") txdb=loadDb(extpath, packageName='GenomicFeatures') IDDB = org.Mm.eg.db CDSdbraw=REFCDS(txdb,IDDB) ## ----eval=TRUE------------------------------------------------------------- DFGENEraw=GENEXP_CDS(CDSdbraw, flsall, Strandtype="forward") ## ----eval=FALSE------------------------------------------------------------ # flsall = dir(bamdir,".bam") # flsall=paste0(bamdir,flsall) # names(flsall)=dir(bamdir,".bam") ## ----sessionInfo----------------------------------------------------------- sessionInfo()