## ----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') ## ----eval=FALSE--------------------------------------------------------------- # library(APAlyzer) ## ---- echo = FALSE------------------------------------------------------------ options(warn=-1) suppressMessages(library(APAlyzer)) ## ----eval=TRUE---------------------------------------------------------------- suppressMessages(library("TBX20BamSubset")) suppressMessages(library("Rsamtools")) flsall = getBamFileList() flsall ## ----eval=TRUE---------------------------------------------------------------- library(repmis) URL="https://github.com/RJWANGbioinfo/PAS_reference_RData/blob/master/" file="mm9_REF.RData" source_data(paste0(URL,file,"?raw=True")) ## ----eval=TRUE---------------------------------------------------------------- head(refUTRraw,2) ## ----eval=TRUE---------------------------------------------------------------- head(dfIPA,2) ## ----eval=TRUE---------------------------------------------------------------- head(dfLE,2) ## ----eval=TRUE---------------------------------------------------------------- URL="https://github.com/RJWANGbioinfo/PAS_reference_RData/blob/master/" file="hg19_REF.RData" source_data(paste0(URL,file,"?raw=True")) ## ----eval=TRUE---------------------------------------------------------------- refUTRraw=refUTRraw[which(refUTRraw$Chrom=='chr19'),] dfIPAraw=dfIPA[which(dfIPA$Chrom=='chr19'),] dfLEraw=dfLE[which(dfLE$Chrom=='chr19'),] PASREF=REF4PAS(refUTRraw,dfIPAraw,dfLEraw) UTRdbraw=PASREF$UTRdbraw dfIPA=PASREF$dfIPA dfLE=PASREF$dfLE ## ----eval=FALSE--------------------------------------------------------------- # ## build Reference ranges for 3'UTR PASs in mouse # download.file(url='ftp://ftp.ensembl.org/pub/release-99/gtf/mus_musculus/Mus_musculus.GRCm38.99.gtf.gz', # destfile='Mus_musculus.GRCm38.99.gtf.gz') # GTFfile="Mus_musculus.GRCm38.99.gtf.gz" # PASREFraw=PAS2GEF(GTFfile) # refUTRraw=PASREFraw$refUTRraw # dfIPAraw=PASREFraw$dfIPA # dfLEraw=PASREFraw$dfLE # PASREF=REF4PAS(refUTRraw,dfIPAraw,dfLEraw) ## ----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): URL="https://github.com/RJWANGbioinfo/PAS_reference_RData/blob/master/" file="mm9_REF.RData" source_data(paste0(URL,file,"?raw=True")) ## ----eval=TRUE---------------------------------------------------------------- #human(hg19): URL="https://github.com/RJWANGbioinfo/PAS_reference_RData/blob/master/" file="hg19_REF.RData" source_data(paste0(URL,file,"?raw=True")) ## ----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) ## ----eval=FALSE--------------------------------------------------------------- # APAVolcano(test_3UTRsing, PAS='3UTR', Pcol = "pvalue", top=5, main='3UTR APA') ## ----out.width = '75%', echo = FALSE------------------------------------------ library(knitr) include_graphics("REDvoca.png") ## ----eval=FALSE--------------------------------------------------------------- # APABox(test_3UTRsing, xlab = "APAreg", ylab = "RED", plot_title = NULL) ## ----out.width = '75%', echo = FALSE------------------------------------------ library(knitr) include_graphics("REDbox.png") ## ---- 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--------------------------------------------------------------- # ## Extract 3 prime most alignment of a paired-end # ## bam file and saved into a new bam file # Bamfile='/path/to/inputdir/input.bam' # Outdir='/path/to/outdir/' # StrandType="forward-reverse" ## "forward-reverse", or "reverse-forward" or "NONE" # ThreeMostPairBam (BamfilePath=Bamfile, # OutDirPath=Outdir, # StrandType='forward-reverse') ## ----eval=FALSE--------------------------------------------------------------- # ThreemostBamfile="test.3most.bam" # IPA_OUTraw=PASEXP_IPA(dfIPA, dfLE, ThreemostBamfile, Strandtype="forward", SeqType="ThreeMostPairEnd") ## ----eval=FALSE--------------------------------------------------------------- # download_testbam() # flsall <- dir(getwd(),".bam") # flsall<-paste0(getwd(),'/',flsall) # names(flsall)<-gsub('.bam','',dir(getwd(),".bam")) ## ----eval=FALSE--------------------------------------------------------------- # library(repmis) # URL="https://github.com/RJWANGbioinfo/PAS_reference_RData/blob/master/" # file="mm9_REF.RData" # source_data(paste0(URL,file,"?raw=True")) # PASREF=REF4PAS(refUTRraw,dfIPA,dfLE) # UTRdbraw=PASREF$UTRdbraw # dfIPA=PASREF$dfIPA # dfLE=PASREF$dfLE ## ----eval=FALSE--------------------------------------------------------------- # UTR_APA_OUT=PASEXP_3UTR(UTRdbraw, flsall, Strandtype="invert") # IPA_OUT=PASEXP_IPA(dfIPA, dfLE, flsall, Strandtype="invert", nts=1) ## ----eval=FALSE--------------------------------------------------------------- # ############# 3utr APA ################# # sampleTable = data.frame(samplename = c('Heart_rep1', # 'Heart_rep2', # 'Heart_rep3', # 'Heart_rep4', # 'Testis_rep1', # 'Testis_rep2', # 'Testis_rep3', # 'Testis_rep4'), # condition = c(rep("Heart",4), # rep("Testis",4))) # # test_3UTRAPA=APAdiff(sampleTable,UTR_APA_OUT, # conKET='Heart', # trtKEY='Testis', # PAS='3UTR', # CUTreads=5) # ## ---- echo = FALSE------------------------------------------------------------ options(warn=-1) extpath = system.file("extdata", "mm9_tissues_rawout.RData", package="APAlyzer") load(extpath) ## ----eval=TRUE---------------------------------------------------------------- table(test_3UTRAPA$APAreg) APAVolcano(test_3UTRAPA, PAS='3UTR', Pcol = "pvalue", plot_title='3UTR APA') APABox(test_3UTRAPA, xlab = "APAreg", ylab = "RED", plot_title = NULL) ## ----eval=FALSE--------------------------------------------------------------- # ############# IPA ################# # test_IPA=APAdiff(sampleTable, # IPA_OUT, # conKET='Heart', # trtKEY='Testis', # PAS='IPA', # CUTreads=5) ## ----eval=TRUE---------------------------------------------------------------- table(test_IPA$APAreg) APAVolcano(test_IPA, PAS='IPA', Pcol = "pvalue", plot_title='IPA') APABox(test_IPA, xlab = "APAreg", ylab = "RED", plot_title = NULL) ## ----eval=FALSE--------------------------------------------------------------- # flsall = dir(bamdir,".bam") # flsall=paste0(bamdir,flsall) # names(flsall)=dir(bamdir,".bam") ## ----sessionInfo-------------------------------------------------------------- sessionInfo()