inst/doc/ASpli.R

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

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


###################################################
### code chunk number 2: loadAspli
###################################################
library(ASpli)
#devtools::load_all("~/Documents/ASpli-labo/")


###################################################
### code chunk number 3: quickStart1 (eval = FALSE)
###################################################
## genome   <- loadDb("txdb.sqlite")
## features <- binGenome(genome)
## targets  <- data.frame(bam = c("CT_1.BAM", "CT_2.BAM","CT_3.BAM",
##                                "TR_1.BAM", "TR_2.BAM", "TR_3.BAM"),
##                   genotype = c( "CT", "CT", "CT",  "TR", "TR", "TR" ),
##                   stringsAsFactors = FALSE )
## mBAMs <- data.frame(bam=c("CT.BAM", "TR.BAM"),condition=c("CT","TR"))


###################################################
### code chunk number 4: quickStart2 (eval = FALSE)
###################################################
## gbcounts  <- gbCounts( features = features, 
##                            targets = targets, 
##                            minReadLength = 100, maxISize = 50000,
##                            libType="SE", 
##                            strandMode=0)


###################################################
### code chunk number 5: quickStart2b (eval = FALSE)
###################################################
## asd         <-  jCounts(counts = gbcounts, 
##                      features = features, 
##                      minReadLength = 125L,
##                      libType="SE", 
##                      strandMode=0)


###################################################
### code chunk number 6: quickStart3 (eval = FALSE)
###################################################
## gb        <- gbDUreport(counts, contrast = c(-1, 1))
## jdur      <- jDUreport(    asd, contrast = c(-1, 1))


###################################################
### code chunk number 7: quickStart4 (eval = FALSE)
###################################################
## sr      <- splicingReport(gb, jdur, counts)
## is      <- integrateSignals(sr,asd)


###################################################
### code chunk number 8: quickStart5 (eval = FALSE)
###################################################
## exportSplicingReports( sr, output.dir="results")
## exportIntegratedSignals(is,sr=sr,
##                           output.dir = "results",
##                           counts=counts,features=features,asd=asd,
##                           mergedBams = mBAMs)


###################################################
### code chunk number 9: installation (eval = FALSE)
###################################################
## if (!requireNamespace("BiocManager", quietly = TRUE))
##     install.packages("BiocManager")
## 
## # The following line initializes usage of Bioc devel branch and should not be necessary after
## # the next official release scheduled for October 2020
## if(Sys.Date()<"2020-11-01") BiocManager::install(version='devel')
## BiocManager::install("ASpli")
## library(ASpli)


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


###################################################
### code chunk number 11: saveTx (eval = FALSE)
###################################################
##  saveDb(TxDB,file="gene.sqlite")


###################################################
### code chunk number 12: binGenome (eval = FALSE)
###################################################
## annFile       <- aspliExampleGTF()
## aTxDb         <- makeTxDbFromGFF(annFile)
## features      <- binGenome( aTxDb ) 
## geneCoord     <- featuresg( features )
## binCoord      <- featuresb( features )
## junctionCoord <- featuresj( features )


###################################################
### code chunk number 13: binGenome2 (eval = FALSE)
###################################################
## symbols       <- data.frame( row.names = genes( aTxDb ), 
##                              symbol = paste( 'This is symbol of gene:',
##                                              genes( aTxDb ) ) )
## features      <- binGenome( aTxDb, geneSymbols = symbols ) 


###################################################
### code chunk number 14: targetsDF2
###################################################
BAMFiles <- c("path_to_bams/CT_time1_rep1.BAM", "path_to_bams/CT_time1_rep2.BAM",
              "path_to_bams/CT_time2_rep1.BAM", "path_to_bams/CT_time2_rep2.BAM",
              "path_to_bams/TR_time1_rep1.BAM", "path_to_bams/TR_time1_rep2.BAM",
              "path_to_bams/TR_time2_rep1.BAM", "path_to_bams/TR_time2_rep2.BAM")
(targets <- data.frame( bam = BAMFiles,
                        genotype = c( 'CT', 'CT', 'CT',  'CT', 
                                      'TR', 'TR', 'TR', 'TR' ),
                        time     = c( 't1', 't1', 't2', 't2', 
                                      't1', 't1', 't2', 't2' ),
                         stringsAsFactors = FALSE ))


###################################################
### code chunk number 15: targetsDF2b
###################################################
getConditions( targets )


###################################################
### code chunk number 16: gbCounts (eval = FALSE)
###################################################
## gbcounts  <- gbCounts( features = features, 
##                            targets = targets, 
##                            minReadLength = 100, maxISize = 50000,
##                            libType="SE", 
##                            strandMode=0)


###################################################
### code chunk number 17: gbCountAccessors (eval = FALSE)
###################################################
## GeneCounts <- countsg(counts)
## GeneRd <- rdsg(counts)
## BinCounts <- countsb(counts)
## BinRd <- rdsb(counts)
## JunctionCounts <- countsj(counts)


###################################################
### code chunk number 18: gbCountsWrite (eval = FALSE)
###################################################
## writeCounts(counts=counts, output.dir = "example")
## writeRds(counts=counts, output.dir = "example")


###################################################
### code chunk number 19: gbCountAccessors2 (eval = FALSE)
###################################################
## e1iCounts <- countse1i(counts)
## ie2Counts <- countsie2(counts)


###################################################
### code chunk number 20: asd (eval = FALSE)
###################################################
## asd         <-  jCounts(counts = gbcounts, 
##                      features = features, 
##                      minReadLength = 100,
##                      libType="SE", 
##                      strandMode=0)


###################################################
### code chunk number 21: asdAccesor (eval = FALSE)
###################################################
## irPIR  <- irPIR( asd )
## altPSI <- altPSI( asd )
## esPSI  <- esPSI( asd )
## allBins      <- joint( asd )
## 
## junctionsPJU <- junctionsPIR( asd )
## junctionsPIR <- junctionsPIR( asd )


###################################################
### code chunk number 22: asdAccesor2 (eval = FALSE)
###################################################
## writeAS(as=asd, output.dir="example")


###################################################
### code chunk number 23: gbDUreport (eval = FALSE)
###################################################
##  gb <- gbDUreport( counts, 
##              minGenReads = 10, 
##              minBinReads = 5,
##              minRds = 0.05, 
##              contrast = NULL, 
##              ignoreExternal = TRUE, 
##              ignoreIo = TRUE, 
##              ignoreI = FALSE,
##              filterWithContrasted = TRUE,
##              verbose = TRUE,
##              formula = NULL,
##              coef = NULL)


###################################################
### code chunk number 24: asdAccesor (eval = FALSE)
###################################################
## geneX  <- genesDE( gb )
## binDU  <- binsDU( gb )


###################################################
### code chunk number 25: gbCountsWrite (eval = FALSE)
###################################################
## writeDU(gb, output.dir = "example")


###################################################
### code chunk number 26: jDUreport (eval = FALSE)
###################################################
## jdu <- jDUreport(asd, 
##             minAvgCounts                       = 5, 
##             contrast                           = NULL,
##             filterWithContrasted               = TRUE,
##             runUniformityTest                  = FALSE,
##             mergedBAMs                         = NULL,
##             maxPValForUniformityCheck          = 0.2,
##             strongFilter                       = TRUE,
##             maxConditionsForDispersionEstimate = 24,
##             formula                            = NULL,
##             coef                               = NULL,
##             maxFDRForParticipation             = 0.05,
##             useSubset                          = FALSE)


###################################################
### code chunk number 27: jDUreportAccesors (eval = FALSE)
###################################################
##  localej( jdu )
##  localec( jdu )
##  
##  anchorj( jdu )
##  anchorc( jdu )
##  
##  jir( jdu )
##  jes( jdu )
##  jalt( jdu )


###################################################
### code chunk number 28: splicingReportA (eval = FALSE)
###################################################
##  binbased( sr )
##  localebased( sr )
##  anchorbased( sr )


###################################################
### code chunk number 29: splicingReportA (eval = FALSE)
###################################################
## writeSplicingReport( sr, output.dir = "test")


###################################################
### code chunk number 30: integrateSignals (eval = FALSE)
###################################################
## is    <- integrateSignals(sr, asd,
##                  bin.FC = 3, bin.fdr = 0.05, bin.inclussion = 0.2,
##                  nonunif = 1, usenonunif = FALSE, 
##                  bjs.inclussion = 10, bjs.fdr = 0.01, 
##                  a.inclussion   = 0.3, a.fdr = 0.01, 
##                  l.inclussion   = 0.3, l.fdr = 0.01,
##                  otherSources = NULL, overlapType = "any")


###################################################
### code chunk number 31: integrateSignalsA (eval = FALSE)
###################################################
##  signals( is )
##  filters( is )


###################################################
### code chunk number 32: exportSplicingReports (eval = FALSE)
###################################################
## exportSplicingReports( sr, 
##                        output.dir="sr",
##                        openInBrowser = FALSE, 
##                        maxBinFDR = 0.2, 
##                        maxJunctionFDR = 0.2 )


###################################################
### code chunk number 33: exportIntegratedSignals (eval = FALSE)
###################################################
## exportIntegratedSignals( is, output.dir="is", 
##                          sr, counts, features, asd,
##                          mergedBams, 
##                          jCompletelyIncluded = FALSE, zoomRegion = 1.5, 
##                          useLog = FALSE, tcex = 1, ntop = NULL, openInBrowser = F, 
##                          makeGraphs = T, bforce=FALSE
##                         )


###################################################
### code chunk number 34: Ex02.a
###################################################
#library( ASpli )
library( GenomicFeatures )
gtfFileName <- aspliExampleGTF()
genomeTxDb <- makeTxDbFromGFF( gtfFileName )
features <- binGenome( genomeTxDb )


###################################################
### code chunk number 35: Ex02.e
###################################################
BAMFiles <- aspliExampleBamList()

targets <- data.frame( 
  row.names = paste0('Sample',c(1:12)),
  bam = BAMFiles,
  f1 = c( 'A','A','A','A','A','A',
          'B','B','B','B','B','B'),
  f2 = c( 'C','C','C','D','D','D',
          'C','C','C','D','D','D'),
  stringsAsFactors = FALSE)


###################################################
### code chunk number 36: Ex02.g
###################################################
getConditions(targets)


###################################################
### code chunk number 37: Ex02.f
###################################################
mBAMs <- data.frame(bam      = sub("_[02]","",targets$bam[c(1,4,7,10)]),
                    condition= c("A_C","A_D","B_C","B_D"))


###################################################
### code chunk number 38: Ex02.i
###################################################
gbcounts  <- gbCounts( features = features, 
                           targets = targets, 
                           minReadLength = 100, maxISize = 50000,
                           libType="SE", 
                           strandMode=0)


###################################################
### code chunk number 39: Ex02.j
###################################################
asd<- jCounts(counts = gbcounts, 
                     features = features, 
                     minReadLength = 100,
                     libType="SE", 
                     strandMode=0,
                      threshold     = 5,
                      minAnchor     = 10)


###################################################
### code chunk number 40: Ex02.k
###################################################
gb      <- gbDUreport(gbcounts,contrast = c( 1, -1, -1, 1 ) )



###################################################
### code chunk number 41: Ex02.k2
###################################################
genesDE(gb)[1:5,]
binsDU(gb)[1:5,]


###################################################
### code chunk number 42: Ex02.l
###################################################
jdur    <- jDUreport(asd, contrast =  c( 1, -1, -1, 1 ))


###################################################
### code chunk number 43: Ex02.kk (eval = FALSE)
###################################################
## localec(jdur)[1:5,]


###################################################
### code chunk number 44: Ex02.kkk
###################################################
localej(jdur)[1:5,1:8]


###################################################
### code chunk number 45: Ex02.m
###################################################
sr      <- splicingReport(gb, jdur, counts = gbcounts)
is      <- integrateSignals(sr,asd)


###################################################
### code chunk number 46: Ex03
###################################################
exportIntegratedSignals( is, output.dir="aspliExample", 
                         sr, gbcounts, features, asd, mBAMs)



###################################################
### code chunk number 47: Ex04
###################################################
 form <- formula(~f1+f2+f1:f2)


###################################################
### code chunk number 48: Ex05
###################################################
 model.matrix(form,targets)


###################################################
### code chunk number 49: Ex02.k (eval = FALSE)
###################################################
## gb      <- gbDUreport(counts, formula = form , coef = 4)
## jdur    <- jDUreport(asd, formula = form, coef = 4 ,
##                      runUniformityTest = TRUE,
##                      mergedBams = mBAMs)


###################################################
### code chunk number 50: Ex02.kp (eval = FALSE)
###################################################
## targetPaired <- targets[c(1,4,7,10),]
## 
## gbcounts  <- gbCounts( features = features, 
##                            targets = targets, 
##                            minReadLength = 100, maxISize = 50000,
##                            libType="SE", 
##                            strandMode=0)
## asd    <- jCounts(counts = gbcounts, 
##                      features = features, 
##                      minReadLength = 100,
##                      libType="SE", 
##                      strandMode=0)
## form   <- formula(~f1+f2)
## gb     <- gbDUreport(gbcounts, formula = form)
## #jdur   <- jDUreport(asd    , formula = form)
## #sr <- splicingReport(gb, jdur, counts = counts)
## #is <- integrateSignals(sr,asd,bjs.fdr = 0.1, bjs.inclussion = 0.1, l.inclussion=0.001,l.fdr = 1)
## # exportSplicingReports(sr,output.dir="paired")
## #exportIntegratedSignals(is,output.dir="paired",sr,counts,features,asd,mBAMs,tcex=2)

Try the ASpli package in your browser

Any scripts or data that you put into this service are public.

ASpli documentation built on Nov. 8, 2020, 5:21 p.m.