###################################################
### chunk number 1: 
###################################################
library(rnaSeqMap)


###################################################
### chunk number 2:  eval=FALSE
###################################################
## xmap.connect("human-24")


###################################################
### chunk number 3: 
###################################################
data(sample_data_rnaSeqMap)
ls()
class(rs.list)
length(rs.list)
names(rs.list)


###################################################
### chunk number 4:  eval=FALSE
###################################################
## rs <- newSeqReads(g.ch,st,en,str)
## rs <- addExperimentsToReadset(rs,1:6)


###################################################
### chunk number 5:  eval=FALSE
###################################################
## rs <- newSeqReadsFromGene(g)
## rs <- addExperimentsToReadset(rs,1:6)


###################################################
### chunk number 6: 
###################################################
rs <- rs.list[[1]]  # first gene from sample data
nd.cov <- getCoverageFromRS(rs,1:6)
distribs(nd.cov)[52000:52010,]    # getting a fragment of the gene


###################################################
### chunk number 7:  eval=FALSE
###################################################
## nsums <- getSumsExp()


###################################################
### chunk number 8: 
###################################################
nsums <- c(1067676,25033137, 85415567,21236253,18152006,25132926)


###################################################
### chunk number 9: 
###################################################
nsums
nd.cov <- normalizeBySum(nd.cov, nsums)
distribs(nd.cov)[52000:52010,]


###################################################
### chunk number 10: 
###################################################
nd.reg <- findRegionsAsND(nd.cov,6,10)
distribs(nd.reg)[52000:52010,]


###################################################
### chunk number 11: 
###################################################
findRegionsAsIR(nd.cov,6,10,1)


###################################################
### chunk number 12: 
###################################################
apply(distribs(nd.cov),2,max)
apply(distribs(nd.cov),2,mean)
apply(distribs(nd.cov),2,sum)
apply(distribs(nd.cov),2,sd)


###################################################
### chunk number 13: 
###################################################
nd.si <- getSIFromND(nd.cov, c(2,3))
distribs(nd.si)[52000:52010,]


###################################################
### chunk number 14:  eval=FALSE
###################################################
## distrSIPlot(nd.cov, 1, 2, mi=5 ,minsup=5)


###################################################
### chunk number 15:  eval=FALSE
###################################################
## plotSI (5,54459000,54461000,-1, 2, 3)


###################################################
### chunk number 16: 
###################################################
nd.fc <- getFCFromND(nd.cov, c(2,3))
distribs(nd.fc)[52000:52010,]


###################################################
### chunk number 17: 
###################################################
nd.rbcov <- regionBasedCoverage(nd.cov, seqq=1:10, maxexp=40, minsup=10)
distribs(nd.rbcov)[52000:52010,]


###################################################
### chunk number 18:  eval=FALSE
###################################################
## distrCOVPlot(nd.cov, c(3,5))


###################################################
### chunk number 19:  eval=FALSE
###################################################
## distrCOVPlot(nd.rbcov, c(3,5))


###################################################
### chunk number 20:  eval=FALSE
###################################################
## my.genes <- geneInChromosome(22, 20000000, 20400000, 1)
## my.genes
## my.spaces <- spaceInChromosome(22, 20000000, 20400000, 1)
## my.spaces


###################################################
### chunk number 21:  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] ))
## }


###################################################
### chunk number 22:  eval=FALSE
###################################################
## plotCoverageHistogram(1,11000,13500,1,1,skip=50)


###################################################
### chunk number 23:  eval=FALSE
###################################################
## distrCOVPlotg(names(rs.list)[16], 1:3)