################################################### ### chunk number 1: options ################################################### #line 41 "vignettes/GenomicRanges/inst/doc/GenomicRangesUseCases.Rnw" options(width=72) ################################################### ### chunk number 2: initialize ################################################### #line 50 "vignettes/GenomicRanges/inst/doc/GenomicRangesUseCases.Rnw" library(GenomicRanges) ################################################### ### chunk number 3: EatonEtAlChIPseq ################################################### #line 66 "vignettes/GenomicRanges/inst/doc/GenomicRangesUseCases.Rnw" library(EatonEtAlChIPseq) data(orcAlignsRep1) orcAlignsRep1 ################################################### ### chunk number 4: EatonAlignmentFiltering-1 ################################################### #line 101 "vignettes/GenomicRanges/inst/doc/GenomicRangesUseCases.Rnw" subsetRep1 <- orcAlignsRep1[alignData(orcAlignsRep1)[["nMismatchBestHit"]] <= 3L] length(subsetRep1) / length(orcAlignsRep1) subsetRep1 <- subsetRep1[occurrenceFilter(withSread=FALSE)(subsetRep1)] length(subsetRep1) / length(orcAlignsRep1) subsetRep1 ################################################### ### chunk number 5: EatonOrcRanges ################################################### #line 118 "vignettes/GenomicRanges/inst/doc/GenomicRangesUseCases.Rnw" rangesRep1 <- as(subsetRep1, "GRanges") head(rangesRep1, 3) ################################################### ### chunk number 6: EatonAddSeqlengths ################################################### #line 126 "vignettes/GenomicRanges/inst/doc/GenomicRangesUseCases.Rnw" seqlengths(rangesRep1) <- 784333 ################################################### ### chunk number 7: EatonNegStarts ################################################### #line 141 "vignettes/GenomicRanges/inst/doc/GenomicRangesUseCases.Rnw" negRangesRep1 <- rangesRep1[strand(rangesRep1) == "-"] negStartsRep1 <- resize(negRangesRep1, 1) ################################################### ### chunk number 8: EatonPossibleEnds ################################################### #line 149 "vignettes/GenomicRanges/inst/doc/GenomicRangesUseCases.Rnw" posRangesRep1 <- rangesRep1[strand(rangesRep1) == "+"] posEndsRep1 <- shift(posRangesRep1, 99) posEndsRep1 <- resize(posEndsRep1, 100) strand(posEndsRep1) <- "-" ################################################### ### chunk number 9: EatonMatchingPairs ################################################### #line 161 "vignettes/GenomicRanges/inst/doc/GenomicRangesUseCases.Rnw" strandMatching <- findOverlaps(negStartsRep1, posEndsRep1) posKeep <- unique(subjectHits(strandMatching)) negKeep <- unique(queryHits(strandMatching)) length(posKeep) / length(posEndsRep1) length(negKeep) / length(negStartsRep1) (length(posKeep) + length(negKeep)) / length(orcAlignsRep1) ################################################### ### chunk number 10: EatonAlignmentFiltering-2 ################################################### #line 170 "vignettes/GenomicRanges/inst/doc/GenomicRangesUseCases.Rnw" posFilteredRangesRep1 <- posRangesRep1[posKeep] negFilteredRangesRep1 <- negRangesRep1[negKeep] ################################################### ### chunk number 11: CoverageWeights ################################################### #line 185 "vignettes/GenomicRanges/inst/doc/GenomicRangesUseCases.Rnw" posWeights <- c(seq(0.01, 1, length = 100), rep(c(1, 0), c(101, 200))) negWeights <- rev(posWeights) plot(-200:200, posWeights, xlab = "Relative Position", ylab = "Coverage Weight", type = "l") ################################################### ### chunk number 12: EatonStartCoverage ################################################### #line 206 "vignettes/GenomicRanges/inst/doc/GenomicRangesUseCases.Rnw" posStartsCoverRep1 <- coverage(resize(posFilteredRangesRep1, 1)) negStartsCoverRep1 <- coverage(resize(negFilteredRangesRep1, 1)) ################################################### ### chunk number 13: EatonExtendedCoverage ################################################### #line 216 "vignettes/GenomicRanges/inst/doc/GenomicRangesUseCases.Rnw" posExtCoverRep1 <- round(runwtsum(posStartsCoverRep1, k = 401, wt = posWeights, endrule = "constant")) negExtCoverRep1 <- round(runwtsum(negStartsCoverRep1, k = 401, wt = negWeights, endrule = "constant")) ################################################### ### chunk number 14: plotCoverage ################################################### #line 231 "vignettes/GenomicRanges/inst/doc/GenomicRangesUseCases.Rnw" plotCoverage <- function(x, xlab = "Position", ylab = "Coverage",...) { plot(c(start(x), length(x)), c(runValue(x), tail(runValue(x), 1)), type = "s", col = "blue", xlab = xlab, ylab = ylab, ...) } plotStrandedCoverage <- function(positive, negative, xlab = "Position", ylab = "Coverage",...) { ylim <- min(max(positive), max(negative)) * c(-1, 1) plotCoverage(positive, ylim = ylim, ...) lines(c(start(negative), length(negative)), - c(runValue(negative), tail(runValue(negative), 1)), type = "s", col = "red") abline(h = 0, col = "dimgray") } ################################################### ### chunk number 15: EatonPlotStrands ################################################### #line 255 "vignettes/GenomicRanges/inst/doc/GenomicRangesUseCases.Rnw" plotStrandedCoverage(posExtCoverRep1[[1]], negExtCoverRep1[[1]]) ################################################### ### chunk number 16: EatonExtendedCoverage ################################################### #line 272 "vignettes/GenomicRanges/inst/doc/GenomicRangesUseCases.Rnw" combExtCoverRep1 <- pmin(posExtCoverRep1, negExtCoverRep1) quantile(combExtCoverRep1, c(0.5, 0.9, 0.95)) ################################################### ### chunk number 17: EatonPeaks ################################################### #line 285 "vignettes/GenomicRanges/inst/doc/GenomicRangesUseCases.Rnw" peaksRep1 <- slice(combExtCoverRep1, lower = 5) peakMaxsRep1 <- viewMaxs(peaksRep1) tail(sort(peakMaxsRep1[[1]]), 30) peaksRep1 <- peaksRep1[peakMaxsRep1 >= 28] peakRangesRep1 <- GRanges("chrXIV", as(peaksRep1[[1]], "IRanges"), seqlengths = seqlengths(rangesRep1)) length(peakRangesRep1) ################################################### ### chunk number 18: EatonExtendedReads ################################################### #line 301 "vignettes/GenomicRanges/inst/doc/GenomicRangesUseCases.Rnw" data(orcPeaksRep1) countOverlaps(orcPeaksRep1, peakRangesRep1) countOverlaps(peakRangesRep1, orcPeaksRep1) ################################################### ### chunk number 19: YeastData ################################################### #line 337 "vignettes/GenomicRanges/inst/doc/GenomicRangesUseCases.Rnw" library(Rsamtools) testFile <- system.file("bam", "isowt5_13e.bam", package="leeBamViews") aligns <- readBamGappedAlignments(testFile) ################################################### ### chunk number 20: GetAnnoations ################################################### #line 363 "vignettes/GenomicRanges/inst/doc/GenomicRangesUseCases.Rnw" library(GenomicFeatures) txdb <- makeTranscriptDbFromUCSC(genome="sacCer2", tablename="sgdGene") ################################################### ### chunk number 21: exonsBy ################################################### #line 373 "vignettes/GenomicRanges/inst/doc/GenomicRangesUseCases.Rnw" exonRanges <- exonsBy(txdb, "tx") length(exonRanges) exonRanges[1] ################################################### ### chunk number 22: ammendData ################################################### #line 387 "vignettes/GenomicRanges/inst/doc/GenomicRangesUseCases.Rnw" seqlevels(exonRanges) newlvls <- c(paste("chr", as.roman(1:16), sep=""), "chrM", "2micron") newlvls seqlevels(exonRanges) <- newlvls # reorder the levels names(newlvls) <- seqlevels(aligns) newlvls seqlevels(aligns) <- newlvls # rename the levels seqlengths(aligns) <- seqlengths(exonRanges) ################################################### ### chunk number 23: count ################################################### #line 407 "vignettes/GenomicRanges/inst/doc/GenomicRangesUseCases.Rnw" counts <- countOverlaps(exonRanges, aligns) ################################################### ### chunk number 24: numBases ################################################### #line 418 "vignettes/GenomicRanges/inst/doc/GenomicRangesUseCases.Rnw" numBases <- sum(width(exonRanges)) geneLengthsInKB <- numBases / 1000 ################################################### ### chunk number 25: millionsMapped ################################################### #line 429 "vignettes/GenomicRanges/inst/doc/GenomicRangesUseCases.Rnw" millionsMapped <- sum(counts) / 1000000 ################################################### ### chunk number 26: RPM ################################################### #line 435 "vignettes/GenomicRanges/inst/doc/GenomicRangesUseCases.Rnw" # counted reads / total reads in millions rpm <- counts / millionsMapped ################################################### ### chunk number 27: RPKM ################################################### #line 442 "vignettes/GenomicRanges/inst/doc/GenomicRangesUseCases.Rnw" # reads per million per geneLength in Kb rpkm <- rpm / geneLengthsInKB ################################################### ### chunk number 28: DESeqFrame eval=FALSE ################################################### ## #line 457 "vignettes/GenomicRanges/inst/doc/GenomicRangesUseCases.Rnw" ## deframe <- data.frame(counts, counts2) ################################################### ### chunk number 29: sort ################################################### #line 475 "vignettes/GenomicRanges/inst/doc/GenomicRangesUseCases.Rnw" sortedRPKM <- sort(rpkm) highScoreGenes <- tail(sortedRPKM) ################################################### ### chunk number 30: annotate1 ################################################### #line 487 "vignettes/GenomicRanges/inst/doc/GenomicRangesUseCases.Rnw" txs <- transcripts(txdb, vals=list(tx_id=names(highScoreGenes)), columns=c("tx_id","gene_id")) systNames <- as.vector(unlist(elementMetadata(txs)["gene_id"])) ################################################### ### chunk number 31: annnotate2 ################################################### #line 498 "vignettes/GenomicRanges/inst/doc/GenomicRangesUseCases.Rnw" library(org.Sc.sgd.db) toTable(org.Sc.sgdGENENAME[systNames]) ################################################### ### chunk number 32: filterKnowns ################################################### #line 522 "vignettes/GenomicRanges/inst/doc/GenomicRangesUseCases.Rnw" filtData <- subsetByOverlaps(aligns, exonRanges) length(filtData) ################################################### ### chunk number 33: filterKnowns2 ################################################### #line 538 "vignettes/GenomicRanges/inst/doc/GenomicRangesUseCases.Rnw" filtData2 <- subsetByOverlaps(aligns, transcriptsBy(txdb, "gene")) length(filtData2) ################################################### ### chunk number 34: coverage ################################################### #line 557 "vignettes/GenomicRanges/inst/doc/GenomicRangesUseCases.Rnw" cov <- coverage(filtData) ################################################### ### chunk number 35: coverageSubset ################################################### #line 568 "vignettes/GenomicRanges/inst/doc/GenomicRangesUseCases.Rnw" cov <- cov[13] ################################################### ### chunk number 36: islands ################################################### #line 575 "vignettes/GenomicRanges/inst/doc/GenomicRangesUseCases.Rnw" islands <- slice(cov, lower = 1) ################################################### ### chunk number 37: continous >1000 ################################################### #line 582 "vignettes/GenomicRanges/inst/doc/GenomicRangesUseCases.Rnw" transcribedRegions <- islands[width(islands) > 1000] txr <- islands[width(islands) > 1000] ################################################### ### chunk number 38: getSequence ################################################### #line 593 "vignettes/GenomicRanges/inst/doc/GenomicRangesUseCases.Rnw" library(BSgenome.Scerevisiae.UCSC.sacCer2) getYeastSequence <- function(data) { chr <- rep("chrXIII",length(start(data[[1]]))) starts <- start(data[[1]]) ends <- end(data[[1]]) strands <- rep("+",length(start(data[[1]]))) getSeq(Scerevisiae, names=chr, start=starts, end=ends, strand=strands) } DNASet <- DNAStringSet(getYeastSequence(transcribedRegions)) ################################################### ### chunk number 39: SessionInfo ################################################### #line 638 "vignettes/GenomicRanges/inst/doc/GenomicRangesUseCases.Rnw" sessionInfo()