## ----eval=FALSE------------------------------------------------------------
#  	#ALIGN:
#  		sample.fastq -> sample.sam
#  	#CONVERT:
#  		samtools view -b -S sample.sam -o sample.bam
#  	#SORT:
#  		samtools sort sample.bam sample_sorted
#  	#INDEXING:
#  		samtools index sample_sorted.bam

## ----eval=TRUE-------------------------------------------------------------
library(wavClusteR)
filename <- system.file( "extdata", "example.bam", package = "wavClusteR" )
Bam <- readSortedBam(filename = filename)
Bam

## ----eval=TRUE-------------------------------------------------------------
countTable <- getAllSub( Bam, minCov = 10 )
head( countTable )

## ----fig.width=5, fig.height=5, fig.align='center', eval=TRUE--------------
plotSubstitutions( countTable, highlight = "TC" )

## ----eval=FALSE------------------------------------------------------------
#  model <- fitMixtureModel(countTable, substitution = "TC")

## ----eval=TRUE-------------------------------------------------------------
data(model)
str(model)

## ----fig.width=7, fig.height=4.5, fig.align='center', eval=TRUE------------
(support <- getExpInterval( model, bayes = TRUE ) )

## ----fig.width=7, fig.height=4.5, fig.align='center', eval=TRUE------------
(support <- getExpInterval( model, bayes = FALSE, leftProb = 0.9, rightProb = 	0.9 ) )

## ----fig.width=7, fig.height=5, fig.align='center', eval=TRUE--------------
plotSubstitutions( countTable, highlight = "TC", model )

## ----eval=TRUE-------------------------------------------------------------
highConfSub <- getHighConfSub( countTable, 
                               support = support, 
                               substitution = "TC" )
head( highConfSub )                               

## ----eval=FALSE------------------------------------------------------------
#  highConfSub <- getHighConfSub( countTable,
#                                 supportStart = 0.2,
#                                 supportEnd = 0.7,
#                                 substitution = "TC" )
#  head( highConfSub )

## ----eval=TRUE-------------------------------------------------------------
coverage <- coverage( Bam )
coverage$chrX

## ----eval=TRUE-------------------------------------------------------------
clusters <- getClusters( highConfSub = highConfSub,
                         		coverage = coverage,
                         		sortedBam = Bam,
                         		method = "mrn",
                         		threshold = 1,
                         		cores = 1 )
clusters

## ----eval=TRUE-------------------------------------------------------------
clusters <- getClusters( highConfSub = highConfSub,
                         		coverage = coverage,
                         		sortedBam = Bam,
                         		method = "mrn",
                         		cores = 1 )
clusters

## ----eval=TRUE-------------------------------------------------------------
require(BSgenome.Hsapiens.UCSC.hg19)

wavclusters <- filterClusters( clusters = clusters, 
                   highConfSub = highConfSub,
                   coverage = coverage, 
                   model = model, 
                   genome = Hsapiens, 
                   refBase = "T", 
                   minWidth = 12)

wavclusters

## ----eval=FALSE------------------------------------------------------------
#  exportHighConfSub( highConfSub = highConfSub,
#                     filename = "hcTC.bed",
#                     trackname = "hcTC",
#                     description = "hcTC" )

## ----eval=FALSE------------------------------------------------------------
#  exportClusters( clusters = wavclusters,
#                  filename = "wavClusters.bed",
#                  trackname = "wavClusters",
#                  description = "wavClusters" )

## ----eval=FALSE------------------------------------------------------------
#  exportCoverage( coverage = coverage, filename = "coverage.bigWig" )

## ----eval=FALSE------------------------------------------------------------
#  txDB <- makeTxDbFromUCSC(genome = "hg19", tablename = "ensGene")

## ----eval=FALSE------------------------------------------------------------
#  annotateClusters( clusters = wavclusters,
#                txDB = txDB,
#                plot = TRUE,
#                verbose = TRUE)

## ----eval=FALSE------------------------------------------------------------
#  getMetaGene( clusters = wavclusters,
#               txDB = txDB,
#               upstream = 1e3,
#               downstream = 1e3,
#               nBins = 40,
#               nBinsUD = 10,
#               minLength = 1,
#               plot = TRUE,
#               verbose = TRUE )

## ----eval=FALSE------------------------------------------------------------
#  getMetaTSS( sortedBam = Bam,
#              txDB = txDB,
#              upstream = 1e3,
#              downstream = 1e3,
#              nBins = 40,
#              unique = FALSE,
#              plot = TRUE,
#              verbose = TRUE )

## ----fig.width=5, fig.height=5, fig.align='center', eval=TRUE--------------
plotSizeDistribution( clusters = wavclusters, showCov = TRUE, col = "skyblue2" )

## ----fig.width=5, fig.height=5, fig.align='center', eval=FALSE-------------
#  plotStatistics( clusters = wavclusters,
#                  corMethod = "spearman",
#                  lower = panel.smooth )

## ----eval=TRUE-------------------------------------------------------------
sessionInfo()