### R code from vignette source 'ASpli.Rnw'

###################################################
### code chunk number 1: style
###################################################
BiocStyle::latex()
options(continue=" ")


###################################################
### code chunk number 2: installation (eval = FALSE)
###################################################
## source("https://bioconductor.org/biocLite.R")
## biocLite("ASpli")


###################################################
### code chunk number 3: loadASpli
###################################################
library(ASpli)


###################################################
### code chunk number 4: makeTx (eval = FALSE)
###################################################
## TxDb <- makeTxDbFromGFF(
##   file="genes.gtf",
##   format="gtf")


###################################################
### code chunk number 5: binGenome (eval = FALSE)
###################################################
## features <- binGenome(TxDb) 
## GenesCoord <- featuresg(features)
## BinsCoord <- featuresb(features)
## JunctionsCoord <- featuresj(features)


###################################################
### code chunk number 6: targetsDF
###################################################
bamFiles <- c("Ct1.bam", "Ct2.bam", "Mut1.bam","Mut2.bam" )
targets <- data.frame(bam=bamFiles,
                    condition=rep(c("CT", "Mut"), each=2), 
                    row.names=sub("\\.bam$", "", bamFiles))


###################################################
### code chunk number 7: targetsDF2
###################################################
bkg <- rep(c("CT", "Mut"), each=4)
time <- rep(rep(c("T1", "T2"), each=2),2)
bamFiles <- c("Ct_time1_rep1.bam", "Ct_time1_rep2.bam",
         "Ct_time2_rep1.bam", "Ct_time2_rep2.bam",
         "Mut_time1_rep1.bam","Mut_time1_rep2.bam",
         "Mut_time2_rep1.bam","Mut_time2_rep2.bam")
targets_2 <- data.frame(bam=bamFiles,
                        condition=paste(bkg, time, sep="."),
                        row.names=sub("\\.bam$", "", bamFiles))


###################################################
### code chunk number 8: loadBam (eval = FALSE)
###################################################
## bam <- loadBAM(targets)


###################################################
### code chunk number 9: loadBam (eval = FALSE)
###################################################
## counts <- readCounts(features, bam, l=100L, targets = targets, maxISize=5000)


###################################################
### code chunk number 10: writeFamil1 (eval = FALSE)
###################################################
## GeneCounts <- countsg(counts)
## GeneRd <- rdsg(counts)
## 
## BinCounts <- countsb(counts)
## BinRd <- rdsb(counts)
## 
## JunctionCounts <- countsj(counts)
## e1iCounts <- countse1i(counts)
## ie2Counts <- countsie2(counts)
## 
## writeCounts(counts=counts, output.dir = "example")
## writeRds(counts=counts, output.dir = "example")


###################################################
### code chunk number 11: AsDiscover (eval = FALSE)
###################################################
## pair <- c("CT","KD")
## as <- AsDiscover(counts, targets, features, bam, l=100L, pair=pair)


###################################################
### code chunk number 12: writeAS (eval = FALSE)
###################################################
## irPIR <- irPIR(as)
## altPSI <- altPSI(as)
## esPSI <- esPSI(as)
## junctionsPIR <- junctionsPIR(as)
## junctionsPSI <- junctionsPSI(as)
## 
## writeAS(as=as, output.dir="example")


###################################################
### code chunk number 13: group (eval = FALSE)
###################################################
## pair <- c("CT","KD")
## group <- c(rep("CT", 4),rep("KD", 4))
## du <- DUreport(counts, targets, pair, group)


###################################################
### code chunk number 14: DEXSeq (eval = FALSE)
###################################################
##  du_DEXSeq <- DUreport_DEXSeq(counts, targets, pair, group)


###################################################
### code chunk number 15: export (eval = FALSE)
###################################################
## writeDU(du, output.dir="example")
## genesde <- genesDE(du)
## binsdu <- binsDU(du)
## junctionsdu <- junctionsDU(du)


###################################################
### code chunk number 16: write (eval = FALSE)
###################################################
## writeCounts(counts, "example_counts")
## writeRds(counts, "example_rds")
## writeDU(du, output.dir="example_du");
## writeAS(as=as, output.dir="example_as");
## writeAll(counts=counts, du=du, as=as, output.dir="example_all")


###################################################
### code chunk number 17: targetsPlot
###################################################
bamFiles <- c("CT.bam","KD.bam")
targetsPlot <- data.frame(bam=bamFiles, 
                        sample=sub("\\.bam$", "", bamFiles), 
                        color=c("blue", "black"), 
                        stringsAsFactors=FALSE)


###################################################
### code chunk number 18: plotTopTags (eval = FALSE)
###################################################
## plotTopTags(auxdf, 
##             TxDb, 
##             targetsPlot, 
##             output.dir="testPlots")


###################################################
### code chunk number 19: librariesEx
###################################################
library(RNAseqData.HNRNPC.bam.chr14)


###################################################
### code chunk number 20: loadDb
###################################################
library( AnnotationDbi )
chr14 <- system.file("extdata","chr14.sqlite", package="ASpli")
genome <- loadDb(chr14)


###################################################
### code chunk number 21: binGenome
###################################################
features <- binGenome(genome) 


###################################################
### code chunk number 22: targetsEx
###################################################
targets <- data.frame(bam=RNAseqData.HNRNPC.bam.chr14_BAMFILES,
                       condition=c(rep("CT",4),rep("KD",4)))



###################################################
### code chunk number 23: loadBam
###################################################
bam <- loadBAM(targets)


###################################################
### code chunk number 24: todosJuntos
###################################################
counts <- readCounts(features, bam, l=100L, targets = targets, maxISize=50000)


###################################################
### code chunk number 25: DU
###################################################
pair <- c("CT","KD")
group <- c(rep("CT", 4),rep("KD", 4))
du_HNRNPC <- DUreport(counts, targets, pair, group)  


###################################################
### code chunk number 26: AS
###################################################
as_HNRNPC <- AsDiscover(counts=counts, targets=targets, features=features, bam=bam, l=100L, pair=pair)


###################################################
### code chunk number 27: topTags
###################################################
binsdu <- binsDU(du_HNRNPC)
topTagsBins <- which(binsdu$bin.fdr <= 0.1 & 
                 abs(binsdu$logFC) >=0.58)


###################################################
### code chunk number 28: targetsPlot
###################################################
targetsPlot <- data.frame(bam=targets$bam, 
                        sample=targets$condition, 
                        color=c(rep("blue", 4),rep("red", 4)), 
                        stringsAsFactors=FALSE)


###################################################
### code chunk number 29: auxdf
###################################################
auxdf<-binsdu[topTagsBins,]
#for simplicity, we choose one: LRR1:E005
plotTopTags(auxdf["LRR1:E005",], 
            genome, 
            targetsPlot, 
            output.dir="testPlots")


###################################################
### code chunk number 30: check
###################################################
binsdu <- binsDU(du_HNRNPC)
topTagsBins <- which(binsdu$bin.fdr <= 0.1 & 
                     abs(binsdu$logFC) >=0.58)


###################################################
### code chunk number 31: DEXSeq
###################################################
du_DEXSeq <- DUreport_DEXSeq(counts, targets, pair, group)
binsdx <- binsDU(du_DEXSeq)
topTagsBinsDx <- which(binsdx$bin.fdr <= 0.1 & 
                     abs(binsdx$logFC) >=0.58)


###################################################
### code chunk number 32: venn
###################################################
library(limma)
all <- union(topTagsBins, topTagsBinsDx)
auxdf <- data.frame(edgeR=all%in%topTagsBins,
                    DEXSeq=all%in%topTagsBinsDx)
vennDiagram(auxdf)


###################################################
### code chunk number 33: sessionInfo
###################################################
sessionInfo()