################################################### ### chunk number 1: b1 ################################################### #line 114 "vignettes/BSgenome/inst/doc/GenomeSearching.Rnw" library(BSgenome.Celegans.UCSC.ce2) ls("package:BSgenome.Celegans.UCSC.ce2") Celegans ################################################### ### chunk number 2: b2 ################################################### #line 121 "vignettes/BSgenome/inst/doc/GenomeSearching.Rnw" class(Celegans) ################################################### ### chunk number 3: b3 ################################################### #line 131 "vignettes/BSgenome/inst/doc/GenomeSearching.Rnw" organism(Celegans) provider(Celegans) providerVersion(Celegans) seqnames(Celegans) mseqnames(Celegans) ################################################### ### chunk number 4: b4 ################################################### #line 143 "vignettes/BSgenome/inst/doc/GenomeSearching.Rnw" Celegans$chrI ################################################### ### chunk number 5: b5 ################################################### #line 175 "vignettes/BSgenome/inst/doc/GenomeSearching.Rnw" chrI <- Celegans$chrI length(chrI) ################################################### ### chunk number 6: b6 ################################################### #line 181 "vignettes/BSgenome/inst/doc/GenomeSearching.Rnw" afI <- alphabetFrequency(chrI) afI sum(afI) == length(chrI) ################################################### ### chunk number 7: b7 ################################################### #line 188 "vignettes/BSgenome/inst/doc/GenomeSearching.Rnw" p1 <- "ACCCAGGGC" countPattern(p1, chrI) ################################################### ### chunk number 8: b8 ################################################### #line 197 "vignettes/BSgenome/inst/doc/GenomeSearching.Rnw" countPattern(p1, chrI, max.mismatch=1) ################################################### ### chunk number 9: b9 ################################################### #line 203 "vignettes/BSgenome/inst/doc/GenomeSearching.Rnw" m1 <- matchPattern(p1, chrI, max.mismatch=1) m1[4:6] class(m1) ################################################### ### chunk number 10: b10 ################################################### #line 211 "vignettes/BSgenome/inst/doc/GenomeSearching.Rnw" mismatch(p1, m1[4:6]) ################################################### ### chunk number 11: b11 ################################################### #line 220 "vignettes/BSgenome/inst/doc/GenomeSearching.Rnw" p2 <- DNAString("AAGCCTAAGCCTAAGCCTAA") m2 <- matchPattern(p2, chrI, max.mismatch=2) m2[1:4] p2 == m2[1:4] mismatch(p2, m2[1:4]) ################################################### ### chunk number 12: b12 ################################################### #line 230 "vignettes/BSgenome/inst/doc/GenomeSearching.Rnw" m2[p2 == m2] m2[p2 != m2] ################################################### ### chunk number 13: c1 ################################################### #line 301 "vignettes/BSgenome/inst/doc/GenomeSearching.Rnw" ce2dict0_file <- system.file("extdata", "ce2dict0.fa", package="BSgenome") ce2dict0 <- read.DNAStringSet(ce2dict0_file, "fasta") ce2dict0 ################################################### ### chunk number 14: c2 ################################################### #line 308 "vignettes/BSgenome/inst/doc/GenomeSearching.Rnw" writeHits <- function(seqname, matches, strand, file="", append=FALSE) { if (file.exists(file) && !append) warning("existing file ", file, " will be overwritten with 'append=FALSE'") if (!file.exists(file) && append) warning("new file ", file, " will have no header with 'append=TRUE'") hits <- data.frame(seqname=rep.int(seqname, length(matches)), start=start(matches), end=end(matches), strand=rep.int(strand, length(matches)), patternID=names(matches), check.names=FALSE) write.table(hits, file=file, append=append, quote=FALSE, sep="\t", row.names=FALSE, col.names=!append) } runAnalysis1 <- function(dict0, outfile="") { library(BSgenome.Celegans.UCSC.ce2) seqnames <- seqnames(Celegans) seqnames_in1string <- paste(seqnames, collapse=", ") cat("Target:", providerVersion(Celegans), "chromosomes", seqnames_in1string, "\n") append <- FALSE for (seqname in seqnames) { subject <- Celegans[[seqname]] cat(">>> Finding all hits in chromosome", seqname, "...\n") for (i in seq_len(length(dict0))) { patternID <- names(dict0)[i] pattern <- dict0[[i]] plus_matches <- matchPattern(pattern, subject) names(plus_matches) <- rep.int(patternID, length(plus_matches)) writeHits(seqname, plus_matches, "+", file=outfile, append=append) append <- TRUE rcpattern <- reverseComplement(pattern) minus_matches <- matchPattern(rcpattern, subject) names(minus_matches) <- rep.int(patternID, length(minus_matches)) writeHits(seqname, minus_matches, "-", file=outfile, append=append) } cat(">>> DONE\n") } } ################################################### ### chunk number 15: c3 ################################################### #line 382 "vignettes/BSgenome/inst/doc/GenomeSearching.Rnw" runAnalysis1(ce2dict0, outfile="ce2dict0_ana1.txt") ################################################### ### chunk number 16: c4 ################################################### #line 390 "vignettes/BSgenome/inst/doc/GenomeSearching.Rnw" hits1 <- read.table("ce2dict0_ana1.txt", header=TRUE) nrow(hits1) ################################################### ### chunk number 17: c5 ################################################### #line 396 "vignettes/BSgenome/inst/doc/GenomeSearching.Rnw" table(hits1$seqname) ################################################### ### chunk number 18: c6 ################################################### #line 401 "vignettes/BSgenome/inst/doc/GenomeSearching.Rnw" hits1_table <- table(hits1$patternID) hits1_table ################################################### ### chunk number 19: c7 ################################################### #line 407 "vignettes/BSgenome/inst/doc/GenomeSearching.Rnw" hits1_table[hits1_table == max(hits1_table)] # pattern(s) with more hits ################################################### ### chunk number 20: c8 ################################################### #line 412 "vignettes/BSgenome/inst/doc/GenomeSearching.Rnw" setdiff(names(ce2dict0), hits1$patternID) # pattern(s) with no hits ################################################### ### chunk number 21: c9 ################################################### #line 417 "vignettes/BSgenome/inst/doc/GenomeSearching.Rnw" plotGenomeHits <- function(bsgenome, seqnames, hits) { chrlengths <- seqlengths(bsgenome)[seqnames] XMAX <- max(chrlengths) YMAX <- length(seqnames) plot.new() plot.window(c(1, XMAX), c(0, YMAX)) axis(1) axis(2, at=seq_len(length(seqnames)), labels=rev(seqnames), tick=FALSE, las=1) ## Plot the chromosomes for (i in seq_len(length(seqnames))) lines(c(1, chrlengths[i]), c(YMAX + 1 - i, YMAX + 1 - i), type="l") ## Plot the hits for (i in seq_len(nrow(hits))) { seqname <- hits$seqname[i] y0 <- YMAX + 1 - match(seqname, seqnames) if (hits$strand[i] == "+") { y <- y0 + 0.05 col <- "red" } else { y <- y0 - 0.05 col <- "blue" } lines(c(hits$start[i], hits$end[i]), c(y, y), type="l", col=col, lwd=3) } } ################################################### ### chunk number 22: c10 eval=FALSE ################################################### ## #line 446 "vignettes/BSgenome/inst/doc/GenomeSearching.Rnw" ## plotGenomeHits(Celegans, seqnames(Celegans), hits1) ################################################### ### chunk number 23: d1 ################################################### #line 469 "vignettes/BSgenome/inst/doc/GenomeSearching.Rnw" library(hgu95av2probe) tmpseq <- DNAStringSet(hgu95av2probe$sequence) someStats <- function(v) { GC <- DNAString("GC") CG <- DNAString("CG") sapply(seq_len(length(v)), function(i) { y <- v[[i]] c(alphabetFrequency(y)[1:4], GC=countPattern(GC, y), CG=countPattern(CG, y)) } ) } someStats(tmpseq[1:10]) ################################################### ### chunk number 24: f1 ################################################### #line 530 "vignettes/BSgenome/inst/doc/GenomeSearching.Rnw" library(BSgenome.Hsapiens.UCSC.hg19) chrY <- Hsapiens$chrY chrY chrM <- Hsapiens$chrM chrM ################################################### ### chunk number 25: f2 ################################################### #line 582 "vignettes/BSgenome/inst/doc/GenomeSearching.Rnw" active(masks(chrY))["RM"] <- TRUE chrY ################################################### ### chunk number 26: f3 ################################################### #line 595 "vignettes/BSgenome/inst/doc/GenomeSearching.Rnw" active(masks(chrY)) <- FALSE active(masks(chrY))["AGAPS"] <- TRUE alphabetFrequency(unmasked(chrY)) alphabetFrequency(chrY) ################################################### ### chunk number 27: f4 ################################################### #line 610 "vignettes/BSgenome/inst/doc/GenomeSearching.Rnw" as(chrY, "XStringViews") ################################################### ### chunk number 28: f5 ################################################### #line 616 "vignettes/BSgenome/inst/doc/GenomeSearching.Rnw" gaps(as(chrY, "XStringViews")) ################################################### ### chunk number 29: f6 ################################################### #line 621 "vignettes/BSgenome/inst/doc/GenomeSearching.Rnw" width(gaps(as(chrY, "XStringViews"))) ################################################### ### chunk number 30: f7 ################################################### #line 628 "vignettes/BSgenome/inst/doc/GenomeSearching.Rnw" gaps(chrY) alphabetFrequency(gaps(chrY)) ################################################### ### chunk number 31: f8 ################################################### #line 635 "vignettes/BSgenome/inst/doc/GenomeSearching.Rnw" af0 <- alphabetFrequency(unmasked(chrY)) af1 <- alphabetFrequency(chrY) af2 <- alphabetFrequency(gaps(chrY)) all(af0 == af1 + af2) ################################################### ### chunk number 32: f9 ################################################### #line 643 "vignettes/BSgenome/inst/doc/GenomeSearching.Rnw" active(masks(chrY)) <- TRUE af1 <- alphabetFrequency(chrY) af1 gaps(chrY) af2 <- alphabetFrequency(gaps(chrY)) af2 all(af0 == af1 + af2) ################################################### ### chunk number 33: f10 ################################################### #line 661 "vignettes/BSgenome/inst/doc/GenomeSearching.Rnw" Ebox <- "CANNTG" active(masks(chrY)) <- FALSE countPattern(Ebox, chrY, fixed=FALSE) ################################################### ### chunk number 34: f11 ################################################### #line 673 "vignettes/BSgenome/inst/doc/GenomeSearching.Rnw" countPattern(Ebox, chrY, fixed=c(pattern=FALSE,subject=TRUE)) ################################################### ### chunk number 35: f12 ################################################### #line 680 "vignettes/BSgenome/inst/doc/GenomeSearching.Rnw" active(masks(chrY))[c("AGAPS", "AMB")] <- TRUE alphabetFrequency(chrY, baseOnly=TRUE) # no ambiguities countPattern(Ebox, chrY, fixed=FALSE) ################################################### ### chunk number 36: f13 ################################################### #line 687 "vignettes/BSgenome/inst/doc/GenomeSearching.Rnw" chr2 <- Hsapiens$chr2 active(masks(chr2))[-2] <- FALSE alphabetFrequency(gaps(chr2)) ################################################### ### chunk number 37: e1 ################################################### #line 734 "vignettes/BSgenome/inst/doc/GenomeSearching.Rnw" runOneStrandAnalysis <- function(dict0, bsgenome, seqnames, strand, outfile="", append=FALSE) { cat("\nTarget: strand", strand, "of", providerVersion(bsgenome), "chromosomes", paste(seqnames, collapse=", "), "\n") if (strand == "-") dict0 <- reverseComplement(dict0) pdict <- PDict(dict0) for (seqname in seqnames) { subject <- bsgenome[[seqname]] cat(">>> Finding all hits in strand", strand, "of chromosome", seqname, "...\n") mindex <- matchPDict(pdict, subject) matches <- extractAllMatches(subject, mindex) writeHits(seqname, matches, strand, file=outfile, append=append) append <- TRUE cat(">>> DONE\n") } } runAnalysis2 <- function(dict0, outfile="") { library(BSgenome.Celegans.UCSC.ce2) seqnames <- seqnames(Celegans) runOneStrandAnalysis(dict0, Celegans, seqnames, "+", outfile=outfile, append=FALSE) runOneStrandAnalysis(dict0, Celegans, seqnames, "-", outfile=outfile, append=TRUE) } ################################################### ### chunk number 38: e2 ################################################### #line 766 "vignettes/BSgenome/inst/doc/GenomeSearching.Rnw" ce2dict0cw15 <- DNAStringSet(ce2dict0, end=15) ################################################### ### chunk number 39: e3 ################################################### #line 772 "vignettes/BSgenome/inst/doc/GenomeSearching.Rnw" runAnalysis2(ce2dict0cw15, outfile="ce2dict0cw15_ana2.txt") ################################################### ### chunk number 40: g1 ################################################### #line 782 "vignettes/BSgenome/inst/doc/GenomeSearching.Rnw" sessionInfo()