################################################### ### 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: addingFeatureData eval=FALSE ################################################### ## if(require("pd.mapping50k.xba240")){ ## system.time(tmp <- chromosome(sample.snpset)) ## object.size(sample.snpset) ## featureData(sample.snpset)$chromosome <- chromosome(sample.snpset) ## featureData(sample.snpset)$position <- position(sample.snpset) ## system.time(tmp <- chromosome(sample.snpset)) ## } ################################################### ### chunk number 6: 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 7: plot1 ################################################### print(graph.par) ################################################### ### chunk number 8: ################################################### graph.par$cytoband.side <- 3 graph.par$heights <- rev(graph.par$heights) ################################################### ### chunk number 9: plot1a eval=FALSE ################################################### ## ##FIXME ## show(graph.par) ################################################### ### chunk number 10: gw.params ################################################### graph.par <- plot(sample.snpset[, 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 11: plot2 ################################################### print(graph.par) ################################################### ### chunk number 12: graph.par3 ################################################### graph.par <- plot(sample.snpset[chromosome(sample.snpset) %in% c(1, 7, 16, 19, "X"), 2]) graph.par$cex <- 0.8 graph.par$mar <- rep(0.5, 4) graph.par$pch <- c(20, 21, 20) graph.par$bty <- "o" graph.par$cex.axis <- 1.2 graph.par$cex.lab <- 1.5 graph.par$xaxs <- "r" graph.par$abline <- TRUE; graph.par$abline.col <- "black" ################################################### ### chunk number 13: plotSnpChromosomes ################################################### print(graph.par) ################################################### ### chunk number 14: 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 15: plotP ################################################### print(graph.par) ################################################### ### chunk number 16: 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 17: 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 18: cytoband ################################################### plotCytoband("1") ################################################### ### chunk number 19: smoothingExample ################################################### sim1 <- sample.snpset[chromosome(sample.snpset) %in% 1:5, 1] sim1 <- sim1[chromosome(sim1) == "1", ] sim1 <- sim1[order(position(sim1)), ] copyNumber(sim1)[101:150, 1] <- copyNumber(sim1)[101:150, 1] - 1 calls(sim1)[101:150, 1] <- 1 smoothSet <- smoothSnp(sim1, 1:5, 1:3, span=1/10) highlight <- calls(smoothSet)[, 1] <= 0.1 & copyNumber(smoothSet)[, 1] <= 1.5 ################################################### ### chunk number 20: plotSmooth ################################################### op <- par(las=1, mfrow=c(1, 1), mar=c(5, 4, 0.5, 0.5), oma=rep(0, 4)) plot(calls(smoothSet)[, 1], copyNumber(smoothSet)[, 1], ylim=range(copyNumber(smoothSet)), pch=".", cex=3, xlab="% heterozygous calls", ylab="smooth copy number", xaxt="n", xlim=c(-0.05, 30/70+.2)) axis(1, at=pretty(calls(smoothSet)), labels=pretty(calls(smoothSet))) points(calls(smoothSet)[highlight, 1], copyNumber(smoothSet)[highlight, 1], pch=20, col="royalblue", bg="white") par(op) ################################################### ### chunk number 21: plot3 ################################################### graph.par$cex <- 1 graph.par$use.chromosome.size <- TRUE graph.par$main <- "Chromosome 1: Example title" graph.par$label.chromosome <- FALSE print(graph.par) ################################################### ### chunk number 22: summary ################################################### x <- summary(sample.snpset, digits=1) str(x) ################################################### ### chunk number 23: 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 24: chromosomeAnnotation ################################################### list.files(system.file("hg18", package="SNPchip")) ##data(chromosomeAnnotation) ##chromosomeAnnotation[1:5, ] ##data(cytoband) ##cytoband[1:5, ] ################################################### ### chunk number 25: annotationSlot eval=FALSE ################################################### ## annotation(sample.snpset) ## library("pd.mapping50k.xba240") ################################################### ### chunk number 26: getSnpAnnotation eval=FALSE ################################################### ## featureData(sample.snpset) <- getSnpAnnotation(sample.snpset) ## fvarLabels(sample.snpset) ################################################### ### chunk number 27: 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 28: nsnpsInRegion eval=FALSE ################################################### ## library(RSNPper) ## (dbId <- dbSnpId(annSnpset)[snps[2] == featureNames(annSnpset)]) ## dbId <- strsplit(dbId, "rs")[[1]][2] ## print(SNPinfo(dbId)) ################################################### ### chunk number 29: ################################################### toLatex(sessionInfo())