################################################### ### 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())