################################################### ### chunk number 1: ################################################### #line 146 "vignettes/rnaSeqMap/inst/doc/rnaSeqMap.Rnw" library(rnaSeqMap) ################################################### ### chunk number 2: eval=FALSE ################################################### ## #line 150 "vignettes/rnaSeqMap/inst/doc/rnaSeqMap.Rnw" ## xmap.connect("human-24") ################################################### ### chunk number 3: ################################################### #line 157 "vignettes/rnaSeqMap/inst/doc/rnaSeqMap.Rnw" data(sample_data_rnaSeqMap) ls() class(rs.list) length(rs.list) names(rs.list) ################################################### ### chunk number 4: eval=FALSE ################################################### ## #line 173 "vignettes/rnaSeqMap/inst/doc/rnaSeqMap.Rnw" ## rs <- newSeqReads(g.ch,st,en,str) ## rs <- addExperimentsToReadset(rs,1:6) ################################################### ### chunk number 5: eval=FALSE ################################################### ## #line 180 "vignettes/rnaSeqMap/inst/doc/rnaSeqMap.Rnw" ## rs <- newSeqReadsFromGene(g) ## rs <- addExperimentsToReadset(rs,1:6) ################################################### ### chunk number 6: ################################################### #line 191 "vignettes/rnaSeqMap/inst/doc/rnaSeqMap.Rnw" setSpecies("homo_sapiens") ################################################### ### chunk number 7: ################################################### #line 194 "vignettes/rnaSeqMap/inst/doc/rnaSeqMap.Rnw" rs <- rs.list[[1]] # first gene from sample data rs <- newSeqReads('chr2',220238268,220264744,-1) rs <- newSeqReads(rnaSeqMap:::.chr.convert(2),220143989,220151622,1) #rs <- newSeqReads('chr2',220114433,220142892,-1) f <- c("test1.bam","test2.bam","test3.bam","test4.bam","test5.bam") ff <- sapply(f, function(x) system.file("extdata", x, package="rnaSeqMap")) rs <- getBamData(rs, 1:5, files=ff) nd.cov <- getCoverageFromRS(rs, 1:5) dd <- RleList2matrix(distribs(nd.cov)) dd[990:1000, ] # getting a fragment of the gene ################################################### ### chunk number 8: eval=FALSE ################################################### ## #line 212 "vignettes/rnaSeqMap/inst/doc/rnaSeqMap.Rnw" ## nsums <- getSumsExp() ################################################### ### chunk number 9: ################################################### #line 216 "vignettes/rnaSeqMap/inst/doc/rnaSeqMap.Rnw" nsums <- c(25918,32577,63250,123445,126977) ################################################### ### chunk number 10: ################################################### #line 220 "vignettes/rnaSeqMap/inst/doc/rnaSeqMap.Rnw" nsums nd.cov <- normalizeBySum(nd.cov, nsums) RleList2matrix(distribs(nd.cov))[990:1000,] ################################################### ### chunk number 11: ################################################### #line 246 "vignettes/rnaSeqMap/inst/doc/rnaSeqMap.Rnw" nd.reg <- findRegionsAsND(nd.cov,6,10) RleList2matrix(distribs(nd.reg))[990:1000,] ################################################### ### chunk number 12: ################################################### #line 254 "vignettes/rnaSeqMap/inst/doc/rnaSeqMap.Rnw" findRegionsAsIR(nd.cov,6,10,1) ################################################### ### chunk number 13: ################################################### #line 269 "vignettes/rnaSeqMap/inst/doc/rnaSeqMap.Rnw" dd <- RleList2matrix(distribs(nd.cov)) apply(dd,2,max) apply(dd,2,mean) apply(dd,2,sum) apply(dd,2,sd) ################################################### ### chunk number 14: eval=FALSE ################################################### ## #line 288 "vignettes/rnaSeqMap/inst/doc/rnaSeqMap.Rnw" ## nd.si <- getSIFromND(nd.cov, c(2,3)) ## RleList2matrix(distribs(nd.si))[990:1000,] ## ################################################### ### chunk number 15: ################################################### #line 295 "vignettes/rnaSeqMap/inst/doc/rnaSeqMap.Rnw" distrCOVPlot(nd.cov, 1) ################################################### ### chunk number 16: ################################################### #line 303 "vignettes/rnaSeqMap/inst/doc/rnaSeqMap.Rnw" distrSIPlot(nd.cov, 1,3,10,10) ################################################### ### chunk number 17: eval=FALSE ################################################### ## #line 309 "vignettes/rnaSeqMap/inst/doc/rnaSeqMap.Rnw" ## plotSI(2,223436000,223437000,-1,2,4) ## ################################################### ### chunk number 18: ################################################### #line 319 "vignettes/rnaSeqMap/inst/doc/rnaSeqMap.Rnw" nd.fc <- getFCFromND(nd.cov, c(2,3)) RleList2matrix(distribs(nd.fc))[990:1000,] ################################################### ### chunk number 19: ################################################### #line 331 "vignettes/rnaSeqMap/inst/doc/rnaSeqMap.Rnw" nd.rbcov <- regionBasedCoverage(nd.cov, seqq=1:10, maxexp=40, minsup=10) RleList2matrix(distribs(nd.rbcov))[990:1000,] ################################################### ### chunk number 20: ################################################### #line 340 "vignettes/rnaSeqMap/inst/doc/rnaSeqMap.Rnw" distrCOVPlot(nd.cov, c(1,3)) ################################################### ### chunk number 21: ################################################### #line 343 "vignettes/rnaSeqMap/inst/doc/rnaSeqMap.Rnw" distrCOVPlot(nd.rbcov, c(1,3)) ################################################### ### chunk number 22: eval=FALSE ################################################### ## #line 363 "vignettes/rnaSeqMap/inst/doc/rnaSeqMap.Rnw" ## my.genes <- geneInChromosome(22, 20000000, 20400000, 1) ## my.genes ## my.spaces <- spaceInChromosome(22, 20000000, 20400000, 1) ## my.spaces ################################################### ### chunk number 23: eval=FALSE ################################################### ## #line 372 "vignettes/rnaSeqMap/inst/doc/rnaSeqMap.Rnw" ## test.gene <- function(g, exps, nsums, mi=5, minsup=15) ## { ## rs <- newSeqReadsFromGene(g) ## rs <- addExperimentsToReadset(rs,exps) ## cat("added reads..\n") ## nd.cov <- getCoverageFromRS(rs,exps) ## cat("calculated coverage...\n") ## nd.cov <- normalizeBySum(nd.cov, nsums) ## cat("normalized...\n") ## nd.reg <- findRegionsAsND(nd.cov,as.integer(mi),minsup=minsup) ## ir.reg <- findRegionsAsIR(nd.cov,as.integer(mi),minsup=minsup) ## cat("ran region search algorithm...\n") ## print(ir.reg) ## ## out <- g ## out <- c(out, apply(distribs(nd.cov),2,max)) ## out <- c(out, apply(distribs(nd.cov),2,mean)) ## out <- c(out, apply(distribs(nd.reg),2,max)) ## out ## } ## ## ## test.space <- function(exps, ch,st, en, str, nsums, mi=5, minsup=15) ## { ## g.ch <- rnaSeqMap:::.chromosome.number(ch) ## rs <- newSeqReads(g.ch,st,en,str) ## rs <- addExperimentsToReadset(rs,exps) ## nd.cov <- getCoverageFromRS(rs,exps) ## nd.cov <- normalizeBySum(nd.cov, nsums) ## cat("Running ch ",g.ch," start",st,"\n") ## nd.reg <- findRegionsAsND(nd.cov,as.integer(mi), minsup=minsup) ## ## out <- c(ch,st, en, str) ## out <- c(out, apply(distribs(nd.cov),2,max)) ## out <- c(out, apply(distribs(nd.cov),2,mean)) ## out <- c(out, apply(distribs(nd.reg),2,max)) ## out ## } ## ## interesting.genes <- NULL ## for (i in 1:length(my.genes)) ## { ## cat ("Running gene ", i , "----------------------\n") ## interesting.genes <- rbind(interesting.genes, test.gene(my.genes[i], 1:6, nsums)) ## } ## ## interesting.spaces <- NULL ## for (i in 104:(dim(my.spaces))[1]) ## { ## cat ("Running space ", i , "----------------------\n") ## interesting.spaces <- rbind(interesting.spaces, test.space(1:2, 22, my.spaces[i,1], my.spaces[i,2],my.spaces[i,3] )) ## } ################################################### ### chunk number 24: eval=FALSE ################################################### ## #line 440 "vignettes/rnaSeqMap/inst/doc/rnaSeqMap.Rnw" ## plotCoverageHistogram(2,223435255,223521056,-1,1,500) ## ################################################### ### chunk number 25: eval=FALSE ################################################### ## #line 446 "vignettes/rnaSeqMap/inst/doc/rnaSeqMap.Rnw" ## distrCOVPlotg(names(rs.list)[16], 1:3)