###################################################
### chunk number 1: 
###################################################
library(BSgenome.Hsapiens.UCSC.hg18)


###################################################
### chunk number 2:  eval=FALSE
###################################################
## ## start with 1MB of chr17
## a = BSgenome:::getSeq(Hsapiens,'chr17',start=40000000,width=1000000)
## ## Generate a little random sequence for insertion
## b = paste(sample(c('A','C','T','G'),10000,replace=TRUE),collapse="")
## ## Now, paste together subsets of the sequences to
## ## generate the abnormal chromosome 17 segment
## newseq = paste(substr(a,1,20000),
##   # start deletions
##   substr(a,20020,22000),
##   substr(a,22040,24000),
##   substr(a,24060,26000),
##   substr(a,26080,28000),
##   substr(a,28100,30000),
##   substr(a,30200,32000),
##   substr(a,34000,40000),
##   # an intrachromosomal translocation
##   substr(a,400000,500000),
##   # a large deletion
##   substr(a,70000,100000),
##   # a high-copy number gain
##   substr(a,100001,110000),
##   substr(a,100001,110000),
##   substr(a,100001,110000),
##   substr(a,100001,110000),
##   substr(a,100001,110000),
##   substr(a,110001,399999),
##   # and a few insertions
##   substr(a,500001,600000),
##   substr(b,1,20),
##   substr(a,600001,602000),
##   substr(b,1,100),
##   substr(a,602001,604000),
##   substr(b,1,1000),
##   substr(a,604001,1000000)
##   )
## f = file("/tmp/sample.sequence.fa","w")
## writeLines('>sample',f)
## writeLines(newseq,f)
## close(f)


###################################################
### chunk number 3: 
###################################################
require(Rsamtools)
scanBamToRangedData <- function(...) {
  tmp = scanBam(...)[[1]]
  toKeep=TRUE
  ranges = IRanges(start=tmp$pos[toKeep],width=nchar(tmp$seq[toKeep]))
  rangeddata = RangedData(ranges,isize=tmp$isize[toKeep],mrnm=tmp$mrnm[toKeep],
    strand=tmp$strand[toKeep],space=tmp$rname[toKeep],flag=tmp$flag[toKeep])
  return(rangeddata)
}


###################################################
### chunk number 4: 
###################################################
fname = system.file('extdata/reads.sorted.bam',package='StructuralVariant')
totalRD = scanBamToRangedData(fname,
  param=ScanBamParam(flag=scanBamFlag(isUnmappedQuery=FALSE)))
totalRD


###################################################
### chunk number 5: 
###################################################
table(space(totalRD))


###################################################
### chunk number 6: 
###################################################
require(bitops)
# find those reads with unmapped mates
# The flag for an unmapped mate is 0x0008
summary(bitAnd(totalRD$flag,0x0008)==0x0008)


###################################################
### chunk number 7: 
###################################################
chr17RD = totalRD['chr17']


###################################################
### chunk number 8: 
###################################################
median(abs(totalRD$isize),na.rm=TRUE)
mad(abs(totalRD$isize),na.rm=TRUE)


###################################################
### chunk number 9: 
###################################################
extendReadsRD <- function(rangeddata,width=100) {
  s = start(rangeddata)
  e = end(rangeddata)
  w = width-width(rangeddata)
  strand = rangeddata$strand
  s[strand=='-']=s[strand=='-']-w[strand=='-']
  e[strand=='+']=e[strand=='+']+w[strand=='+']
  start(rangeddata) <- s
  end(rangeddata) <- e
  return(rangeddata)
}


###################################################
### chunk number 10: 
###################################################
smallInsertRD = chr17RD[(bitAnd(chr17RD$flag,0x0008)!=0x0008) & 
  (chr17RD$mrnm=='chr17') & (abs(chr17RD$isize)<137),]
smallInsertRD = extendReadsRD(smallInsertRD,135)
cvgSI = coverage(smallInsertRD)
sliceSI = slice(cvgSI,10)


###################################################
### chunk number 11: 
###################################################
largeInsertRD = chr17RD[(bitAnd(chr17RD$flag,0x0008)!=0x0008) & 
  (chr17RD$mrnm=='chr17') & (abs(chr17RD$isize)>263),]
largeInsertRD = extendReadsRD(largeInsertRD,135)
cvgLI = coverage(largeInsertRD)
sliceLI = slice(cvgLI,4)


###################################################
### chunk number 12: 
###################################################
unmappedRD = chr17RD[(bitAnd(chr17RD$flag,0x0008)==0x0008),]
unmappedRD = extendReadsRD(unmappedRD,135)
cvgUM = coverage(unmappedRD)
sliceUM = slice(cvgUM,5)


###################################################
### chunk number 13: 
###################################################
sliceLI[['chr17']]


###################################################
### chunk number 14: 
###################################################
plotSV <- function(rangeddata,cvgslice,space='chr17',...) {
  xl = min(start(rangeddata)[space(rangeddata)==space])
  xr = max(end(rangeddata)[space(rangeddata)==space])
  par(las=2)
  plot(c(xl,xr),c(0,3),pch='',xlab="genomic coordinates",ylab="",axes=FALSE,...)
  for(i in 1:length(cvgslice[[space]])) {
    rect(start(cvgslice[[space]])[i],0,end(cvgslice[[space]])[i],0.5,col='red',border='red')
  }
  subrd = rangeddata[rangeddata %in% RangedData(
    IRanges(start=start(cvgslice[[space]]),end=end(cvgslice[[space]])),
    space=rep(space,length(cvgslice[[space]]))),]
  for(i in 1:length(start(subrd))) {
    strandq=ifelse(bitAnd(subrd$flag[i],0x0010)==0x0010,'-','+')
    strandmate=ifelse(bitAnd(subrd$flag[i],0x0020)==0x0020,'-','+')
    if(strandq=='+' & strandmate=='-') col='red'
    if(strandq=='-' & strandmate=='+') col='green'
    if(strandq=='+' & strandmate=='+') col='blue'
    if(strandq=='-' & strandmate=='-') col='orange'
    lines(list(x=c(start(subrd)[i],start(subrd)[i]+subrd$isize[i]),y=c(0.5,2.5)),col=col)
  }
  axis(side=1)
  box()
}


###################################################
### chunk number 15: 
###################################################
pdf('fullRegionLI.pdf',width=9,height=7)
plotSV(chr17RD,sliceLI,xlim=c(4e7,4.1e7))
dev.off()
pdf('zoomDeletionsLI.pdf',width=9,height=7)
plotSV(chr17RD,sliceLI,xlim=c(4.002e7,4.0045e7))
dev.off()


###################################################
### chunk number 16: 
###################################################
toLatex(sessionInfo())