Bioconductor package for analyzing 4C-seq data

### 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(GenomicAlignments)
libraryFile <- system.file("extdata", "vfl_aagctt_catg_mm9_54_vp.csv", package="Basic4Cseq")
bamFile <- system.file("extdata", "fetalLiverShort.bam", package="Basic4Cseq")
liverReads <- readGAlignments(bamFile)
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()