contiBAIT.R
In contiBAIT: Improves Early Build Genome Assemblies using Strand-Seq Data

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

###################################################
### code chunk number 1: strandSeqFreqTableExamplea
###################################################
# Read in BAM files. Path denotes location of the BAM files.
# Returns a vector of file locations

library(contiBAIT)
bamFileList <- list.files(
path=file.path(system.file(package='contiBAIT'), 'extdata'),
pattern=".bam$",
full.names=TRUE)


###################################################
### code chunk number 2: makeChrTableExample
###################################################
# build chr table from BAM file in bamFileList

exampleChrTable <- makeChrTable(bamFileList[1]) 

exampleChrTable


###################################################
### code chunk number 3: makeChrTableExampleb
###################################################

exampleDividedChr <- makeChrTable(bamFileList[1], splitBy=1000000)
exampleDividedChr


###################################################
### code chunk number 4: mapGapOverlapExampleb
###################################################

library(rtracklayer)
# Download GRCh38/hg38 gap track from UCSC
gapFile <- import.bed("http://genome.ucsc.edu/cgi-bin/hgTables?hgsid=465319523_SLOtFPExny48YZFaXBh4sSTzuMcA&boolshad.hgta_printCustomTrackHeaders=0&hgta_ctName=tb_gap&hgta_ctDesc=table+browser+query+on+gap&hgta_ctVis=pack&hgta_ctUrl=&fbQual=whole&fbUpBases=200&fbDownBases=200&hgta_doGetBed=get+BED")
# Create fake SCE file with four regions that overlap a gap
sceFile <- GRanges(rep('chr4',4), 
IRanges(c(1410000, 1415000, 1420000, 1425000), 
c(1430000, 1435000, 1430000, 1435000)))
overlappingFragments <- mapGapFromOverlap(sceFile,
 gapFile, 
 exampleChrTable, 
 overlapNum=4)
show(overlappingFragments)


###################################################
### code chunk number 5: strandSeqFreqTableExamplea
###################################################
# Create a strandFreqTable instance 

strandFrequencyList <- strandSeqFreqTable(bamFileList, 
filter=exampleDividedChr,
qual=10, 
pairedEnd=FALSE,
BAITtables=TRUE)


###################################################
### code chunk number 6: strandSeqFreqTableExampleb
###################################################
# Returned list consisting of two data.frames
strandFrequencyList

# Exclude frequencies calculated from
# contigs with less than 10 reads

exampleStrandFreq <- strandFrequencyList[[1]]
exampleReadCounts <- strandFrequencyList[[2]]
exampleStrandFreq[which(exampleReadCounts < 10)] <- NA 


###################################################
### code chunk number 7: strandSeqFreqTableExamplec
###################################################

# Assess the quality of the libraries being analysed
plotWCdistribution(exampleStrandFreq)


###################################################
### code chunk number 8: preprocessStrandTableExamplea
###################################################
# Convert strand frequencies to strand calls.

exampleStrandStateMatrix <- preprocessStrandTable(
exampleStrandFreq, 
lowQualThreshold=0.8)

exampleStrandStateMatrix[[1]]


###################################################
### code chunk number 9: clusterContigsExamplea
###################################################
exampleWCMatrix <- exampleStrandStateMatrix[[1]]

clusteredContigs <- clusterContigs(exampleWCMatrix, randomise=FALSE)

reorientedMatrix <- reorientAndMergeLGs(clusteredContigs,
 exampleWCMatrix)

exampleLGList <- reorientedMatrix[[3]]

exampleLGList
exampleLGList[[1]]


###################################################
### code chunk number 10: clusterContigsExampleb
###################################################
plotLGDistances(exampleLGList, exampleWCMatrix)


###################################################
### code chunk number 11: clustercontigsExamplec
###################################################

plotLGDistances(exampleLGList, exampleWCMatrix, lg=1)


###################################################
### code chunk number 12: orderAllLinkageGroupsExample
###################################################
contigOrder <- orderAllLinkageGroups(exampleLGList,
exampleWCMatrix,
exampleStrandFreq,
exampleReadCounts,
whichLG=1,
saveOrdered=TRUE)

contigOrder[[1]]


###################################################
### code chunk number 13: orderAllLinkageGroupsExampleb
###################################################

plotContigOrder(contigOrder[[1]])


###################################################
### code chunk number 14: orderAllLinkageGroupsExamplec
###################################################
contigOrderAll <- orderAllLinkageGroups(exampleLGList,
exampleWCMatrix,
exampleStrandFreq,
exampleReadCounts)

contigOrderAll[[1]]


###################################################
### code chunk number 15: ideogramExample
###################################################
# extract elements from strandSeqFreqTable list
WatsonFreqList <- strandFrequencyList[[3]]
CrickFreqList <- strandFrequencyList[[4]]

# subset elements to only analyze one library
singleWatsonLibrary <- StrandReadMatrix(WatsonFreqList[,2, drop=FALSE])
singleCrickLibrary <- StrandReadMatrix(CrickFreqList[,2, drop=FALSE]) 

# Run ideogram plotter
ideogramPlot(singleWatsonLibrary,
singleCrickLibrary,
exampleDividedChr)


###################################################
### code chunk number 16: ideogramExampleb
###################################################

ideogramPlot(singleWatsonLibrary,
singleCrickLibrary,
exampleDividedChr,
orderFrame=contigOrderAll[[1]])


###################################################
### code chunk number 17: ideogramExamplec
###################################################
ideogramPlot(WatsonFreqList,
CrickFreqList,
exampleDividedChr,
orderFrame=contigOrder[[1]],
plotBy='chr',
showPage=1)


###################################################
### code chunk number 18: writeBedExample
###################################################

writeBed(exampleDividedChr,
reorientedMatrix[[2]],
contigOrder[[1]])


###################################################
### code chunk number 19: makeChrTableExamplec
###################################################
makeBoxPlot(exampleDividedChr, exampleLGList)


###################################################
### code chunk number 20: makeChrTableExampled
###################################################
barplotLinkageGroupCalls(exampleLGList, exampleDividedChr)


###################################################
### code chunk number 21: fig1plot
###################################################
library(diagram)

par(mar = c(1, 1, 1, 1))
 openplotmat()
# elpos <- coordinates (c(5, 5, 5, 5, 5, 5, 5, 5))
 elpos <- coordinates (rep(4,7))

fromto <- matrix(ncol = 2, byrow = TRUE, data = c(
	    3, 6,
 	    3, 8,
 	    6, 5,
 	    6, 7,
 	    6,9,
 	    6,10,
 	    8,7,
 	    8,11,
 	    8,20,
 	    10,14,
 	    14,18,
 	    18,19,
 	    10,11,
 	    10,15,
 	    19,15,
 	    19,20,
 	    19,23,
 	    23,27,
 	    9,17,
 	    14,13,
 	    17,22,
 	    22,23,
 	    18,23,
 	    23,24


))

 nr <- nrow(fromto)
 arrpos <- matrix(ncol = 2, nrow = nr)
 for (i in 1:nr)
arrpos[i, ] <- straightarrow (to = elpos[fromto[i, 2], ], from = elpos[fromto[i, 1], ], lwd = 2, arr.pos = 0.6, arr.width=0.15, arr.length = 0.5)

textround(elpos[3,], 0.04, lab = "BAMFILE", box.col = "grey70", shadow.col = "grey10", shadow.size = 0.005, cex = 0.8)
textrect(elpos[8,], 0.09, 0.04,lab = "make\nChrTable", box.col = "white", shadow.col = "grey10", shadow.size = 0.005, cex = 0.7)
textrect(elpos[6,], 0.09, 0.04,lab = "strandSeq\nFreqTable", box.col = "white", shadow.col = "grey10", shadow.size = 0.005, cex = 0.7)
textrect(elpos[10,], 0.09, 0.04,lab = "preprocess\nStrandTable", box.col = "white", shadow.col = "grey10", shadow.size = 0.005, cex = 0.7)
textrect(elpos[14,], 0.09, 0.04,lab = "cluster\nContigs", box.col = "white", shadow.col = "grey10", shadow.size = 0.005, cex = 0.7)
textrect(elpos[18,], 0.09, 0.04,lab = "reorient\nLinkageGroups", box.col = "white", shadow.col = "grey10", shadow.size = 0.005, cex = 0.7)
textrect(elpos[19,], 0.09, 0.04,lab = "merge\nLinkageGroup", box.col = "white", shadow.col = "grey10", shadow.size = 0.005, cex = 0.7)
textrect(elpos[23,], 0.09, 0.04,lab = "orderAll\nLinkageGroups", box.col = "white", shadow.col = "grey10", shadow.size = 0.005, cex = 0.7)
textround(elpos[27,], 0.04, lab = "writeBed", box.col = "grey70", shadow.col = "grey10", shadow.size = 0.005, cex = 0.8)

texthexa(elpos[7,], 0.09, 0.04,lab = "ideogram\nPlot", box.col = "lightblue", shadow.col = "grey10", shadow.size = 0.005, cex = 0.7)
texthexa(elpos[11,], 0.09, 0.04,lab = "make\nBoxPlot", box.col = "lightblue", shadow.col = "grey10", shadow.size = 0.005, cex = 0.7)
texthexa(elpos[5,], 0.09, 0.04,lab = "plotWC\ndistribution", box.col = "lightblue", shadow.col = "grey10", shadow.size = 0.005, cex = 0.7)
texthexa(elpos[15,], 0.09, 0.04,lab = "plotLG\ndistances", box.col = "lightblue", shadow.col = "grey10", shadow.size = 0.005, cex = 0.7)
texthexa(elpos[20,], 0.09, 0.04,lab = "barplot\nLGCalls", box.col = "lightblue", shadow.col = "grey10", shadow.size = 0.005, cex = 0.7)
texthexa(elpos[24,], 0.09, 0.04,lab = "plotContig\nOrder", box.col = "lightblue", shadow.col = "grey10", shadow.size = 0.005, cex = 0.7)




###################################################
### code chunk number 22: fig1
###################################################
library(diagram)

par(mar = c(1, 1, 1, 1))
 openplotmat()
# elpos <- coordinates (c(5, 5, 5, 5, 5, 5, 5, 5))
 elpos <- coordinates (rep(4,7))

fromto <- matrix(ncol = 2, byrow = TRUE, data = c(
	    3, 6,
 	    3, 8,
 	    6, 5,
 	    6, 7,
 	    6,9,
 	    6,10,
 	    8,7,
 	    8,11,
 	    8,20,
 	    10,14,
 	    14,18,
 	    18,19,
 	    10,11,
 	    10,15,
 	    19,15,
 	    19,20,
 	    19,23,
 	    23,27,
 	    9,17,
 	    14,13,
 	    17,22,
 	    22,23,
 	    18,23,
 	    23,24


))

 nr <- nrow(fromto)
 arrpos <- matrix(ncol = 2, nrow = nr)
 for (i in 1:nr)
arrpos[i, ] <- straightarrow (to = elpos[fromto[i, 2], ], from = elpos[fromto[i, 1], ], lwd = 2, arr.pos = 0.6, arr.width=0.15, arr.length = 0.5)

textround(elpos[3,], 0.04, lab = "BAMFILE", box.col = "grey70", shadow.col = "grey10", shadow.size = 0.005, cex = 0.8)
textrect(elpos[8,], 0.09, 0.04,lab = "make\nChrTable", box.col = "white", shadow.col = "grey10", shadow.size = 0.005, cex = 0.7)
textrect(elpos[6,], 0.09, 0.04,lab = "strandSeq\nFreqTable", box.col = "white", shadow.col = "grey10", shadow.size = 0.005, cex = 0.7)
textrect(elpos[10,], 0.09, 0.04,lab = "preprocess\nStrandTable", box.col = "white", shadow.col = "grey10", shadow.size = 0.005, cex = 0.7)
textrect(elpos[14,], 0.09, 0.04,lab = "cluster\nContigs", box.col = "white", shadow.col = "grey10", shadow.size = 0.005, cex = 0.7)
textrect(elpos[18,], 0.09, 0.04,lab = "reorient\nLinkageGroups", box.col = "white", shadow.col = "grey10", shadow.size = 0.005, cex = 0.7)
textrect(elpos[19,], 0.09, 0.04,lab = "merge\nLinkageGroup", box.col = "white", shadow.col = "grey10", shadow.size = 0.005, cex = 0.7)
textrect(elpos[23,], 0.09, 0.04,lab = "orderAll\nLinkageGroups", box.col = "white", shadow.col = "grey10", shadow.size = 0.005, cex = 0.7)
textround(elpos[27,], 0.04, lab = "writeBed", box.col = "grey70", shadow.col = "grey10", shadow.size = 0.005, cex = 0.8)

texthexa(elpos[7,], 0.09, 0.04,lab = "ideogram\nPlot", box.col = "lightblue", shadow.col = "grey10", shadow.size = 0.005, cex = 0.7)
texthexa(elpos[11,], 0.09, 0.04,lab = "make\nBoxPlot", box.col = "lightblue", shadow.col = "grey10", shadow.size = 0.005, cex = 0.7)
texthexa(elpos[5,], 0.09, 0.04,lab = "plotWC\ndistribution", box.col = "lightblue", shadow.col = "grey10", shadow.size = 0.005, cex = 0.7)
texthexa(elpos[15,], 0.09, 0.04,lab = "plotLG\ndistances", box.col = "lightblue", shadow.col = "grey10", shadow.size = 0.005, cex = 0.7)
texthexa(elpos[20,], 0.09, 0.04,lab = "barplot\nLGCalls", box.col = "lightblue", shadow.col = "grey10", shadow.size = 0.005, cex = 0.7)
texthexa(elpos[24,], 0.09, 0.04,lab = "plotContig\nOrder", box.col = "lightblue", shadow.col = "grey10", shadow.size = 0.005, cex = 0.7)