## ----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,
                        p_adjust_methods="fdr")

## ----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,
                        p_adjust_methods="fdr",
                        MultiTest='unpaired t-test')
head(test_3UTRmuti,2)
table(test_3UTRmuti$APAreg)

## ----eval=TRUE----------------------------------------------------------------
test_IPAsing=APAdiff(sampleTable2,
                        IPA_OUTraw, 
                        conKET='NT',
                        trtKEY='KD',
                        PAS='IPA',
                        CUTreads=0,
                        p_adjust_methods="fdr") 
head(test_IPAsing,2)

## ----eval=TRUE----------------------------------------------------------------
test_IPAmuti=APAdiff(sampleTable1,
                        IPA_OUTraw, 
                        conKET='NT',
                        trtKEY='KD',
                        PAS='IPA',
                        CUTreads=0,
                        p_adjust_methods="fdr",
                        MultiTest='unpaired t-test')
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()

