###################################################
### chunk number 1: options
###################################################
options(width=60)


###################################################
### chunk number 2: package
###################################################
library(SNPchip)


###################################################
### chunk number 3: exampledata
###################################################
data(sample.snpset)
sample.snpset


###################################################
### chunk number 4: featureData eval=FALSE
###################################################
## fD <- new("AnnotatedDataFrame", data=data.frame(row.names=featureNames(sample.snpset)),
##           varMetadata=data.frame(labelDescription=character()))
## featureData(sample.snpset) <- fD


###################################################
### chunk number 5: subsetting
###################################################
snpset <- sample.snpset[chromosome(sample.snpset) %in% as.character(1:3), 1:4]
graph.par <- plot(snpset)
class(graph.par)
graph.par$use.chromosome.size <- TRUE
graph.par$pch <- "."
graph.par$cex <- 3
graph.par$oma <- c(3, 3, 3, 0.5)


###################################################
### chunk number 6: plot1
###################################################
print(graph.par)


###################################################
### chunk number 7: 
###################################################
graph.par$cytoband.side <- 3
graph.par$heights <- rev(graph.par$heights)


###################################################
### chunk number 8: plot1a eval=FALSE
###################################################
## ##FIXME
## show(graph.par)


###################################################
### chunk number 9: gw.params
###################################################
graph.par <- plot(sample.snpset[chromosome(sample.snpset) < 23, 1:3], add.cytoband=FALSE)
graph.par$one.ylim <- TRUE
graph.par$mar <- c(0.1, 0.1, 2, 0.1)
graph.par$oma <- c(3, 4, 2, 1)
graph.par$cex <- 2
graph.par$abline <- TRUE
graph.par$cex.lab <- 0.9
graph.par$add.cytoband <- FALSE


###################################################
### chunk number 10: plot2
###################################################
print(graph.par)


###################################################
### chunk number 11: parm
###################################################
parm <- centromere("1")[1]
##data(chromosomeAnnotation)
##parm <- chromosomeAnnotation["1", "centromereStart"]
snpset <- sample.snpset[chromosome(sample.snpset) == "1" & (position(sample.snpset) < parm), 2]
graph.par <- plot(snpset)
graph.par$use.chromosome.size <- FALSE
graph.par$pch <- 21
graph.par$cex <- 1
graph.par$ylim <- c(0.4, 9)
graph.par$cytoband.ycoords <- c(0.5, 0.6)


###################################################
### chunk number 12: plotP
###################################################
print(graph.par)


###################################################
### chunk number 13: addLegend
###################################################
data(sample.snpset)
x <- sample.snpset[chromosome(sample.snpset) == "1", 1]
gp <- plot(x)
gp$legend <- c("AA", "AB", "BB")
gp$legend.col <- unique(gp$col)
gp$legend.bg <- unique(gp$bg)
gp$pch <- 21; gp$cex <- 0.8
gp$label.cytoband <- TRUE
gp$add.centromere <- FALSE
gp$xlab <- ""
gp$legend.bty="o"
gp$ylim[2] <- 9
print(gp)
##plot(gp, x)
##xlim <- c(0,max(position(x)))
##xlim <- range(position(x))
##plotCytoband("1", xlim=xlim, label.cytoband=TRUE)
##xlim <- c(0, max(position(x)))
##plotCytoband("1", xlim=xlim, label.cytoband=TRUE)


###################################################
### chunk number 14:  eval=FALSE
###################################################
## ##gp <- new("ParSnpSet")
## ##gp <- getPar(gp, x)
## gp$pch <- 21; gp$cex <- 0.8
## ##gp$cytoband.side <- 3; 
## ##gp$heights <- rev(gp$heights)
## gp$label.cytoband <- TRUE
## gp$ylim[2] <- 10
## gp$cytoband.ycoords <- c(0.9, 10)
## gp$add.centromere <- FALSE
## gp
## ##plot(gp, x)
## legend("bottomleft", 
##        pch=21, 
##        col=gp$col, 
##        pt.bg=gp$bg, legend=c("AA", "AB", "BB"), bty="n")


###################################################
### chunk number 15: cytoband
###################################################
plotCytoband("1")


###################################################
### chunk number 16: summary
###################################################
x <- summary(sample.snpset, digits=1)
str(x)


###################################################
### chunk number 17: summaryPlot
###################################################
op <- par(mfrow=c(1,1), mar=c(4, 4, 3, 1), las=1)
boxplot(split(copyNumber(sample.snpset[, 1]), chromosome(sample.snpset)), 
        ylab="copy number", main=sampleNames(sample.snpset)[1], log="y")
par(op)


###################################################
### chunk number 18: chromosomeAnnotation
###################################################
list.files(system.file("hg18", package="SNPchip"))
##data(chromosomeAnnotation)
##chromosomeAnnotation[1:5, ]
##data(cytoband)
##cytoband[1:5, ]


###################################################
### chunk number 19: annotationSlot eval=FALSE
###################################################
## annotation(sample.snpset)
## library("pd.mapping50k.xba240")


###################################################
### chunk number 20: getSnpAnnotation eval=FALSE
###################################################
## featureData(sample.snpset) <- getSnpAnnotation(sample.snpset)
## fvarLabels(sample.snpset)


###################################################
### chunk number 21: netAffxAnnotation eval=FALSE
###################################################
## path <- "http://biostat.jhsph.edu/~iruczins/publications/sm/2006.scharpf.bioinfo"
## try(load(url(paste(path, "/mapping/mapping10k.rda", sep=""))))
## colnames(mapping10k$annotation)


###################################################
### chunk number 22: nsnpsInRegion eval=FALSE
###################################################
## library(RSNPper)
## (dbId <- dbSnpId(annSnpset)[snps[2] == featureNames(annSnpset)])
## dbId <- strsplit(dbId, "rs")[[1]][2]
## print(SNPinfo(dbId))


###################################################
### chunk number 23: 
###################################################
toLatex(sessionInfo())