### R code from vignette source 'vignettes/rnaSeqMap/inst/doc/rnaSeqMap.Rnw' ################################################### ### code chunk number 1: rnaSeqMap.Rnw:153-154 ################################################### library(rnaSeqMap) ################################################### ### code chunk number 2: rnaSeqMap.Rnw:157-158 (eval = FALSE) ################################################### ## xmap.connect("human-24") ################################################### ### code chunk number 3: rnaSeqMap.Rnw:164-169 ################################################### data(sample_data_rnaSeqMap) ls() class(rs.list) length(rs.list) names(rs.list) ################################################### ### code chunk number 4: rnaSeqMap.Rnw:180-182 (eval = FALSE) ################################################### ## rs <- newSeqReads(g.ch,st,en,str) ## rs <- addExperimentsToReadset(rs,1:6) ################################################### ### code chunk number 5: rnaSeqMap.Rnw:187-189 (eval = FALSE) ################################################### ## rs <- newSeqReadsFromGene(g) ## rs <- addExperimentsToReadset(rs,1:6) ################################################### ### code chunk number 6: rnaSeqMap.Rnw:198-200 ################################################### setSpecies("homo_sapiens") ################################################### ### code chunk number 7: rnaSeqMap.Rnw:201-213 ################################################### 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 ################################################### ### code chunk number 8: rnaSeqMap.Rnw:219-220 (eval = FALSE) ################################################### ## nsums <- getSumsExp() ################################################### ### code chunk number 9: rnaSeqMap.Rnw:223-224 ################################################### nsums <- c(25918,32577,63250,123445,126977) ################################################### ### code chunk number 10: rnaSeqMap.Rnw:227-230 ################################################### nsums nd.cov <- normalizeBySum(nd.cov, nsums) RleList2matrix(distribs(nd.cov))[990:1000,] ################################################### ### code chunk number 11: rnaSeqMap.Rnw:253-255 ################################################### nd.reg <- findRegionsAsND(nd.cov,6,10) RleList2matrix(distribs(nd.reg))[990:1000,] ################################################### ### code chunk number 12: rnaSeqMap.Rnw:261-262 ################################################### findRegionsAsIR(nd.cov,6,10,1) ################################################### ### code chunk number 13: rnaSeqMap.Rnw:276-281 ################################################### dd <- RleList2matrix(distribs(nd.cov)) apply(dd,2,max) apply(dd,2,mean) apply(dd,2,sum) apply(dd,2,sd) ################################################### ### code chunk number 14: rnaSeqMap.Rnw:295-298 (eval = FALSE) ################################################### ## nd.si <- getSIFromND(nd.cov, c(2,3)) ## RleList2matrix(distribs(nd.si))[990:1000,] ## ################################################### ### code chunk number 15: rnaSeqMap.Rnw:302-303 ################################################### distrCOVPlot(nd.cov, 1) ################################################### ### code chunk number 16: rnaSeqMap.Rnw:310-311 ################################################### distrSIPlot(nd.cov, 1,3,10,10) ################################################### ### code chunk number 17: rnaSeqMap.Rnw:316-318 (eval = FALSE) ################################################### ## plotSI(2,223436000,223437000,-1,2,4) ## ################################################### ### code chunk number 18: rnaSeqMap.Rnw:326-328 ################################################### nd.fc <- getFCFromND(nd.cov, c(2,3)) RleList2matrix(distribs(nd.fc))[990:1000,] ################################################### ### code chunk number 19: rnaSeqMap.Rnw:338-340 ################################################### nd.rbcov <- regionBasedCoverage(nd.cov, seqq=1:10, maxexp=40, minsup=10) RleList2matrix(distribs(nd.rbcov))[990:1000,] ################################################### ### code chunk number 20: rnaSeqMap.Rnw:347-348 ################################################### distrCOVPlot(nd.cov, c(1,3)) ################################################### ### code chunk number 21: rnaSeqMap.Rnw:350-351 ################################################### distrCOVPlot(nd.rbcov, c(1,3)) ################################################### ### code chunk number 22: rnaSeqMap.Rnw:370-374 (eval = FALSE) ################################################### ## my.genes <- geneInChromosome(22, 20000000, 20400000, 1) ## my.genes ## my.spaces <- spaceInChromosome(22, 20000000, 20400000, 1) ## my.spaces ################################################### ### code chunk number 23: rnaSeqMap.Rnw:379-431 (eval = FALSE) ################################################### ## 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] )) ## } ################################################### ### code chunk number 24: rnaSeqMap.Rnw:447-449 (eval = FALSE) ################################################### ## plotCoverageHistogram(2,223435255,223521056,-1,1,500) ## ################################################### ### code chunk number 25: rnaSeqMap.Rnw:453-454 (eval = FALSE) ################################################### ## distrCOVPlotg(names(rs.list)[16], 1:3)