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

###################################################
### code chunk number 1: vignette.Rnw:36-38
###################################################
options(width=90)
options(continue=" ")


###################################################
### code chunk number 2: preliminaries
###################################################
library(Basic4Cseq)


###################################################
### code chunk number 3: createVirtualFragmentLibrary
###################################################
testGenome = DNAString("ATCATGAAGTACTACATGGCACCATGT")
fragmentData = createVirtualFragmentLibrary(chosenGenome = testGenome, firstCutter = "catg", 
	secondCutter = "gtac", readLength = 2,  chromosomeName = "test", libraryName = "")
head(fragmentData)


###################################################
### code chunk number 4: l2 (eval = FALSE)
###################################################
## library(BSgenome.Mmusculus.UCSC.mm9)
## createVirtualFragmentLibrary(chosenGenome = Mmusculus, firstCutter = "aagctt", secondCutter = 
## 	"catg", readLength = 54,  libraryName = "myFullFetalLiverVFL.csv")


###################################################
### code chunk number 5: loadFetalLiverRawData
###################################################
  library(ShortRead)
  libraryFile <- system.file("extdata", "vfl_aagctt_catg_mm9_54_vp.csv", package="Basic4Cseq")
  bamFile <- system.file("extdata", "fetalLiverShort.bam", package="Basic4Cseq")
  liverReads <- readAligned(bamFile, type = "BAM")
  liverReads


###################################################
### code chunk number 6: addPointsOfInterest
###################################################
  pointsOfInterestFile <- system.file("extdata", "fetalLiverVP.bed", package="Basic4Cseq")
  liverPoints<-readPointsOfInterestFile(pointsOfInterestFile)
  liverPoints


###################################################
### code chunk number 7: initializeObject
###################################################
liverData = Data4Cseq(viewpointChromosome = "10", viewpointInterval = c(20879870, 20882209), 
	readLength = 54, pointsOfInterest = liverPoints, rawReads = liverReads)
liverData


###################################################
### code chunk number 8: initializeObject
###################################################
  rawFragments(liverData)<-readsToFragments(liverData, libraryFile)
  liverData


###################################################
### code chunk number 9: initializeObject
###################################################
  # pick near-cis fragments
  nearCisFragments(liverData)<-chooseNearCisFragments(liverData, 
	regionCoordinates = c(20800000, 21100000))
  head(nearCisFragments(liverData))


###################################################
### code chunk number 10: normalizeObject
###################################################
  # normalization of near-cis data
  library("caTools")
  nearCisFragments(liverData)<-normalizeFragmentData(liverData)
  head(nearCisFragments(liverData))


###################################################
### code chunk number 11: qualityControl
###################################################
  getReadDistribution(liverData, useFragEnds = TRUE, outputName = "")


###################################################
### code chunk number 12: visualizeViewpoint (eval = FALSE)
###################################################
##   # This command creates a near-cis plot of the fetal liver's viewpoint data
##   visualizeViewpoint(liverData, plotFileName = "", mainColour = "blue", 
## 	plotTitle = "Fetal Liver Near-Cis Plot", loessSpan = 0.1, maxY = 6000, 
## 	xAxisIntervalLength = 50000, yAxisIntervalLength = 1000)


###################################################
### code chunk number 13: drawHeatmap (eval = FALSE)
###################################################
##   # This command creates a near-cis heatmap plot of the fetal liver data
##   drawHeatmap(liverData, plotFileName = "", xAxisIntervalLength = 50000, bands = 5)


###################################################
### code chunk number 14: exportWig (eval = FALSE)
###################################################
##   printWigFile(liverData, wigFileName = "fetalLiver.wig")


###################################################
### code chunk number 15: transInteractions (eval = FALSE)
###################################################
##   library("RCircos")
##   transInteractions <- system.file("extdata", "transInteractionData.txt", package="Basic4Cseq")
##   ideogramData <- system.file("extdata", "RCircos_GRCm38_ideogram.csv", package="Basic4Cseq")
##   plotTransInteractions(transInteractions, "10", c(20000100, 20001000), ideogramData, 
## 	PlotColor = "blue", expandBands = TRUE, expansionValue = 1000000, plotFileName = "")


###################################################
### code chunk number 16: checkRestrictionEnzyme (eval = FALSE)
###################################################
##   # The demonstration reads are taken from Stadhouders et al's data,
##   # but additional cutter sequences have been added manually for demonstration purposes
##   fetalLiverCutterFile <- system.file("extdata", "fetalLiverCutter.sam", package="Basic4Cseq")
##   checkRestrictionEnzymeSequence("aagctt", fetalLiverCutterFile, "fetalLiverCutter_filtered.sam")


###################################################
### code chunk number 17: extractBED (eval = FALSE)
###################################################
##   printBEDFragmentLibrary(libraryFile, "BEDLibrary_FL_vp.bed")


###################################################
### code chunk number 18: giveEnzymeSequence
###################################################
  enzymeData <- system.file("extdata", "enzymeData.csv", package="Basic4Cseq")
  giveEnzymeSequence(enzymeData, "NlaIII")


###################################################
### code chunk number 19: simulateDigestion
###################################################
  shortTestGenome = "ATCCATGTAGGCTAAGTACACATGTTAAGGTACAGTACAATTGCACGATCAT"
  fragments = simulateDigestion("catg", "gtac", shortTestGenome)
  head(fragments)


###################################################
### code chunk number 20: drawDigestionFragmentHistogram
###################################################
  # This command creates a histogram plot of virtual library fragment length and frequencies
  drawDigestionFragmentHistogram(fragments)


###################################################
### code chunk number 21: prepare4CseqData (eval = FALSE)
###################################################
##     # BWA, samtools and bedtools must be installed
##     # It is assumed that the example data files (from the package) are in the active directory
##     prepare4CseqData("veryShortExample.fastq", "CATG", "veryShortLib.csv", 
## 	referenceGenome = "veryShortReference.fasta")


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