APAlyzer.R
In APAlyzer: A toolkit for APA analysis using RNA-seq data

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