### R code from vignette source 'vignettes/htSeqTools/inst/doc/htSeqTools.Rnw' ################################################### ### code chunk number 1: import1 ################################################### options(width=70) library(htSeqTools) data(htSample) htSample sapply(htSample,nrow) ################################################### ### code chunk number 2: import2 ################################################### Ctrl=as.data.frame(htSample[[1]]) IP1=as.data.frame(htSample[[2]]) head(Ctrl) head(RangedData(Ctrl)) htSample2 <- RangedDataList(Ctrl=RangedData(Ctrl),IP1=RangedData(IP1)) htSample2 head(htSample2[['Ctrl']]) ################################################### ### code chunk number 3: getmds ################################################### cmds1 <- cmds(htSample,k=2) cmds1 ################################################### ### code chunk number 4: figmds ################################################### plot(cmds1) ################################################### ### code chunk number 5: figmds ################################################### plot(cmds1) ################################################### ### code chunk number 6: qcontrol1 ################################################### ssdCoverage(htSample) giniCoverage(htSample,mc.cores=1) ################################################### ### code chunk number 7: figgini1 ################################################### giniCoverage(htSample[['ctrlBatch2']],mc.cores=1,mk.plot=TRUE) ################################################### ### code chunk number 8: figgini2 ################################################### giniCoverage(htSample[['ipBatch2']],mc.cores=1,mk.plot=TRUE) ################################################### ### code chunk number 9: figgini1 ################################################### giniCoverage(htSample[['ctrlBatch2']],mc.cores=1,mk.plot=TRUE) ################################################### ### code chunk number 10: figgini2 ################################################### giniCoverage(htSample[['ipBatch2']],mc.cores=1,mk.plot=TRUE) ################################################### ### code chunk number 11: qcontrol2 ################################################### contamSample <- RangedData(IRanges(rep(1,200),rep(36,200)),space=rep('chr2L',200),strand='+') contamSample <- rbind(htSample[['ctrlBatch1']],contamSample) nrepeats <- tabDuplReads(contamSample) nrepeats ################################################### ### code chunk number 12: qcontrol3 ################################################### q <- which(cumsum(nrepeats/sum(nrepeats))>.999)[1] q fdrest <- fdrEnrichedCounts(nrepeats, use=1:q, components=1) numRepeArtif <- rownames(fdrest[fdrest$fdrEnriched<0.01,])[1] numRepeArtif ################################################### ### code chunk number 13: figrepeats1 ################################################### plot(fdrest$fdrEnriched,type='l',ylab='',xlab='Number of repeats (r)') lines(fdrest$pdfOverall,col=4) lines(fdrest$pdfH0,col=2,lty=2) legend('topright',c('FDR','P(r | no artifact)','P(r)'),lty=c(1,2,1),col=c(1,2,4)) ################################################### ### code chunk number 14: figrepeats1 ################################################### plot(fdrest$fdrEnriched,type='l',ylab='',xlab='Number of repeats (r)') lines(fdrest$pdfOverall,col=4) lines(fdrest$pdfH0,col=2,lty=2) legend('topright',c('FDR','P(r | no artifact)','P(r)'),lty=c(1,2,1),col=c(1,2,4)) ################################################### ### code chunk number 15: figpreprocess1 ################################################### covbefore <- coverage(htSample[['ipBatch2']]) covbefore <- seqselect(covbefore[['chr2L']],295108,297413) plot(as.integer(covbefore),type='l',ylab='Coverage') ################################################### ### code chunk number 16: preprocess1 ################################################### ip2Aligned <- alignPeaks(htSample[['ipBatch2']],strand='strand', npeaks=100) covafter <- coverage(ip2Aligned) covafter <- seqselect(covafter[['chr2L']],295108,297413) covplus <- coverage(htSample[['ipBatch2']][htSample[['ipBatch2']][['strand']]=='+',]) covplus <- seqselect(covplus[['chr2L']],295108,297413) covminus <- coverage(htSample[['ipBatch2']][htSample[['ipBatch2']][['strand']]=='-',]) covminus <- seqselect(covminus[['chr2L']],295108,297413) ################################################### ### code chunk number 17: figpreprocess2 ################################################### plot(as.integer(covafter),type='l',ylab='Coverage') lines(as.integer(covplus),col='blue',lty=2) lines(as.integer(covminus),col='red',lty=2) ################################################### ### code chunk number 18: figpreprocess1 ################################################### covbefore <- coverage(htSample[['ipBatch2']]) covbefore <- seqselect(covbefore[['chr2L']],295108,297413) plot(as.integer(covbefore),type='l',ylab='Coverage') ################################################### ### code chunk number 19: figpreprocess2 ################################################### plot(as.integer(covafter),type='l',ylab='Coverage') lines(as.integer(covplus),col='blue',lty=2) lines(as.integer(covminus),col='red',lty=2) ################################################### ### code chunk number 20: islandCounts ################################################### ip <- rbind(htSample[[2]],htSample[[4]]) ctrl <- rbind(htSample[[1]],htSample[[3]]) pool <- RangedDataList(ip=ip, ctrl=ctrl) counts <- islandCounts(pool,minReads=10) head(counts) ################################################### ### code chunk number 21: enrichedRegions ################################################### mappedreads <- c(ip=nrow(ip),ctrl=nrow(ctrl)) mappedreads regions <- enrichedRegions(sample1=ip,sample2=ctrl,minReads=10,mappedreads=mappedreads,p.adjust.method='BY',pvalFilter=.05,twoTailed=FALSE) head(regions) nrow(regions) ################################################### ### code chunk number 22: enrichedPeaks ################################################### peaks <- enrichedPeaks(regions, sample1=ip, sample2=ctrl, minHeight=100) peaks ################################################### ### code chunk number 23: mergeRegions ################################################### mergeRegions(peaks, maxDist=300, score='height', aggregateFUN='median') ################################################### ### code chunk number 24: analysis2 ################################################### chrLength <- 500000 names(chrLength) <- c('chr2L') chrregions <- enrichedChrRegions(peaks, chrLength=chrLength, windowSize=10^4-1, fdr=0.05, nSims=1) chrregions[['chr2L']] ################################################### ### code chunk number 25: figanalysis2 ################################################### chrregions <- RangedData(ranges=IRanges(start=c(100000,200000),end=c(100100,210000)),space=c('chr2L','chr2R')) plotChrRegions(regions=chrregions, chrLength=c(chr2L=500000,chr2R=500000)) ################################################### ### code chunk number 26: figanalysis2 ################################################### chrregions <- RangedData(ranges=IRanges(start=c(100000,200000),end=c(100100,210000)),space=c('chr2L','chr2R')) plotChrRegions(regions=chrregions, chrLength=c(chr2L=500000,chr2R=500000)) ################################################### ### code chunk number 27: figplots1 ################################################### set.seed(1) peaksanno <- peaks peaksanno[['start_position']] <- start(peaksanno) + runif(nrow(peaksanno),-500,1000) peaksanno[['end_position']] <- peaksanno[['start_position']] + 500 peaksanno[['distance']] <- peaksanno[['start_position']] - start(peaksanno) peaksanno[['strand']] <- sample(c('+','-'),nrow(peaksanno),replace=TRUE) PeakLocation(peaksanno,peakDistance=1000) ################################################### ### code chunk number 28: figplots1 ################################################### set.seed(1) peaksanno <- peaks peaksanno[['start_position']] <- start(peaksanno) + runif(nrow(peaksanno),-500,1000) peaksanno[['end_position']] <- peaksanno[['start_position']] + 500 peaksanno[['distance']] <- peaksanno[['start_position']] - start(peaksanno) peaksanno[['strand']] <- sample(c('+','-'),nrow(peaksanno),replace=TRUE) PeakLocation(peaksanno,peakDistance=1000) ################################################### ### code chunk number 29: regionscov ################################################### cover <- coverage(ip) rcov <- regionsCoverage(space(regions),start(regions),end(regions),cover=cover) names(rcov) rcov[['views']] gridcov <- gridCoverage(rcov) dim(getCover(gridcov)) getViewsInfo(gridcov) ################################################### ### code chunk number 30: figregCoverage ################################################### ylim <- c(0,max(getViewsInfo(gridcov)[['maxCov']])) plot(gridcov, ylim=ylim,lwd=2) for (i in 1:nrow(regions)) lines(getCover(gridcov)[i,], col='gray', lty=2) ################################################### ### code chunk number 31: figregCoverage ################################################### ylim <- c(0,max(getViewsInfo(gridcov)[['maxCov']])) plot(gridcov, ylim=ylim,lwd=2) for (i in 1:nrow(regions)) lines(getCover(gridcov)[i,], col='gray', lty=2)