###################################################
### chunk number 1: lkd
###################################################
#line 98 "vignettes/leeBamViews/inst/doc/leeViews.Rnw"
library(leeBamViews)  # bam files stored in package
bpaths = dir(system.file("bam", package="leeBamViews"), full=TRUE, patt="bam$")
#
# extract genotype and lane information from filenames
#
gt = gsub(".*/", "", bpaths)
gt = gsub("_.*", "", gt)
lane = gsub(".*(.)$", "\\1", gt)
geno = gsub(".$", "", gt)
#
# format the sample-level information appropriately
#
pd = DataFrame(geno=geno, lane=lane, row.names=paste(geno,lane,sep="."))
prd = new("DataFrame")  # protocol data could go here
#
# create the views object, adding some arbitrary experiment-level information
#
bs1 = BamViews(bamPaths=bpaths, bamSamples=pd, 
        bamExperiment=list(annotation="org.Sc.sgd.db"))
bs1
#
# get some sample-level data
#
bamSamples(bs1)$geno


###################################################
### chunk number 2: lkc
###################################################
#line 129 "vignettes/leeBamViews/inst/doc/leeViews.Rnw"
START=c(861250, 863000)
END=c(862750, 864000)
exc = GRanges(seqnames="Scchr13", IRanges(start=START, end=END), strand="+")
bamRanges(bs1) = exc
bs1


###################################################
### chunk number 3: lkcov
###################################################
#line 142 "vignettes/leeBamViews/inst/doc/leeViews.Rnw"
covex =
  RleList(lapply(bamPaths(bs1), function(x)
                 coverage(readGappedAlignments(x))[["Scchr13"]]))
names(covex) = gsub(".bam$", "", basename(bamPaths(bs1)))
head(covex, 3)


###################################################
### chunk number 4: addso
###################################################
#line 152 "vignettes/leeBamViews/inst/doc/leeViews.Rnw"
library(GenomeGraphs)
cov2baseTrack = function(rle, start, end,
   dp = DisplayPars(type="l", lwd=0.5, color="black"),
   countTx=function(x)log10(x+1)) {
 require(GenomeGraphs)
 if (!is(rle, "Rle")) stop("requires instance of Rle")
 dat = runValue(rle)
 loc = cumsum(runLength(rle))
 ok = which(loc >= start & loc <= end)
 makeBaseTrack(base = loc[ok], value=countTx(dat[ok]),
    dp=dp)
}
trs = lapply(covex, function(x) cov2baseTrack(x, START[1], END[1],
   countTx = function(x)pmin(x, 80)))
ac = as.character
names(trs) = paste(ac(bamSamples(bs1)$geno), ac(bamSamples(bs1)$lane), sep="")
library(biomaRt)
mart = useMart("ensembl", "scerevisiae_gene_ensembl")
gr = makeGeneRegion(START, END, chromosome="XIII",
  strand="+", biomart=mart, dp=DisplayPars(plotId=TRUE,
  idRotation=0, idColor="black"))
trs[[length(trs)+1]] = gr
trs[[length(trs)+1]] = makeGenomeAxis()


###################################################
### chunk number 5: lkf
###################################################
#line 177 "vignettes/leeBamViews/inst/doc/leeViews.Rnw"
print( gdPlot( trs, minBase=START[1], maxBase=END[1]) )


###################################################
### chunk number 6: mkps eval=FALSE
###################################################
## #line 182 "vignettes/leeBamViews/inst/doc/leeViews.Rnw"
## plotStrains = function(bs, query, start, end, snames, mart, chr, strand="+") {
##  filtbs = bs[query, ]
##  cov = lapply(filtbs, coverage)
##  covtrs = lapply(cov, function(x) cov2baseTrack(x[[1]], start, end,
##    countTx = function(x) pmin(x,80)))
##  names(covtrs) = snames
##  gr = makeGeneRegion(start, end, chromosome=chr,
##        strand=strand, biomart=mart, dp=DisplayPars(plotId=TRUE,
##        idRotation=0, idColor="black"))
##  grm = makeGeneRegion(start, end, chromosome=chr,
##        strand="-", biomart=mart, dp=DisplayPars(plotId=TRUE,
##        idRotation=0, idColor="black"))
##  covtrs[[length(covtrs)+1]] = gr
##  covtrs[[length(covtrs)+1]] = makeGenomeAxis()
##  covtrs[[length(covtrs)+1]] = grm
##  gdPlot( covtrs, minBase=start, maxBase=end )
## }


###################################################
### chunk number 7: lkda eval=FALSE
###################################################
## #line 216 "vignettes/leeBamViews/inst/doc/leeViews.Rnw"
## data(leeUnn)
## names(leeUnn)
## leeUnn[1:4,1:8]
## table(leeUnn$study)
## l13 = leeUnn[ leeUnn$chr == 13, ]
## l13d = na.omit(l13[ l13$study == "David", ])
## d13r = GRanges(seqnames="Scchr13", IRanges( l13d$start, l13d$end ),
##   strand=ifelse(l13d$strand==1, "+", ifelse(l13d$strand=="0", "*", "-")))
## elementMetadata(d13r)$name = paste("dav13x", 1:length(d13r), sep=".")
## bamRanges(bs1) = d13r
## d13tab = tabulateReads( bs1 )


###################################################
### chunk number 8: makn
###################################################
#line 232 "vignettes/leeBamViews/inst/doc/leeViews.Rnw"
myrn = GRanges(seqnames="Scchr13",
  IRanges(start=seq(861250, 862750, 100), width=100), strand="+")
elementMetadata(myrn)$name = paste("til", 1:length(myrn), sep=".")
bamRanges(bs1) = myrn
tabulateReads(bs1, "+")


###################################################
### chunk number 9: ddee
###################################################
#line 247 "vignettes/leeBamViews/inst/doc/leeViews.Rnw"
library(org.Sc.sgd.db)
library(IRanges)
c13g = get("13", revmap(org.Sc.sgdCHR))  # all genes on chr13
c13loc = unlist(mget(c13g, org.Sc.sgdCHRLOC))  # their 'start' addresses
c13locend = unlist(mget(c13g, org.Sc.sgdCHRLOCEND))
c13locp = c13loc[c13loc>0]     # confine attention to + strand
c13locendp = c13locend[c13locend>0]
ok = !is.na(c13locp) & !is.na(c13locendp)
c13pr = GRanges(seqnames="Scchr13", IRanges(c13locp[ok], c13locendp[ok]),
    strand="+")   # store and clean up names
elementMetadata(c13pr)$name = gsub(".13$", "", names(c13locp[ok]))
c13pr
c13pro = c13pr[ order(ranges(c13pr)), ]


###################################################
### chunk number 10: dolim
###################################################
#line 266 "vignettes/leeBamViews/inst/doc/leeViews.Rnw"
lim = GRanges(seqnames="Scchr13", IRanges(800000,900000), strand="+")
c13prol = c13pro[ which(c13pro %in% lim ), ]


###################################################
### chunk number 11: getm
###################################################
#line 274 "vignettes/leeBamViews/inst/doc/leeViews.Rnw"
bamRanges(bs1) = c13prol
annotab = tabulateReads(bs1, strandmarker="+")


###################################################
### chunk number 12: gettot
###################################################
#line 301 "vignettes/leeBamViews/inst/doc/leeViews.Rnw"
totcnts = totalReadCounts(bs1)


###################################################
### chunk number 13: lkedd
###################################################
#line 309 "vignettes/leeBamViews/inst/doc/leeViews.Rnw"
library(edgeR)
#
# construct an edgeR container for read counts, including 
#   genotype and region (gene) metadata
#
dgel1 = DGEList( counts=t(annotab)[,-c(1,2)], 
   group=factor(bamSamples(bs1)$geno),
   lib.size=totcnts, genes=colnames(annotab))
#
# compute a dispersion factor for the negative binomial model
#
cd = estimateCommonDisp(dgel1)
#
# test for differential expression between two groups
# for each region
#
et12 = exactTest(cd)
#
# display statistics for the comparison
#
tt12 = topTags(et12)
tt12


###################################################
### chunk number 14: lkma
###################################################
#line 340 "vignettes/leeBamViews/inst/doc/leeViews.Rnw"
plotSmear(cd, cex=.8, ylim=c(-10,10))
text(tt12$table$logConc, tt12$table$logFC+.15, as.character(
  tt12$table$genes), cex=.65)


###################################################
### chunk number 15: lks
###################################################
#line 378 "vignettes/leeBamViews/inst/doc/leeViews.Rnw"
sessionInfo()