[BioC] How to use snapCGH with Nimblegen data
Steven McKinney
smckinney at bccrc.ca
Fri Oct 20 02:15:34 CEST 2006
Hi Romain,
Here is a script and some utility functions
that I use to read in and process NimbleGen aCGH
data. I'm attaching the R scripts and have
copied them as clear text at the end, hopefully
one form or other will get through to you.
These scripts work for the data I am currently
processing, but be aware that NimbleGen formats
can and do change, so you should monitor the
inputs and outputs carefully to ensure things
are working for your formats.
I received a lot of help and code examples
from Natalie Thorne (snapCGH developer currently
at Oxford) and Todd Richmond (Manager of Research
Informatics, NimbleGen) which I modified to work
for our situation. Many many thanks to them for
all their good work and help.
Best of luck! Hope these work for you.
Steven McKinney
Statistician
Molecular Oncology and Breast Cancer Program
British Columbia Cancer Research Centre
email: smckinney at bccrc.ca
tel: 604-675-8000 x7561
BCCRC
Molecular Oncology
675 West 10th Ave, Floor 4
Vancouver B.C.
V5Z 1L3
Canada
-----Original Message-----
From: bioconductor-bounces at stat.math.ethz.ch on behalf of Romain Desprat
Sent: Thu 10/19/2006 6:17 AM
To: bioconductor at stat.math.ethz.ch
Subject: [BioC] How to use snapCGH with Nimblegen data
There is this readPositionalInfo function in the snapCGH package tu use
nimblegen data but I don't know how and where in the script I should use it.
Does someone could send me an example of the script using CGH data and
nimblegen format?
Romain
rdesprat at aecom.yu.edu
_______________________________________________
Bioconductor mailing list
Bioconductor at stat.math.ethz.ch
https://stat.ethz.ch/mailman/listinfo/bioconductor
Search the archives: http://news.gmane.org/gmane.science.biology.informatics.conductor
###---------------------------------------------------------------------
## processNimblegenData.R: Example script for reading in and
## processing data from Nimblegen chips.
## This script used to process chips using the XVGA (768 x 1024)
## 1:2 feature density format, (probes at every 2nd location, approx
## 380000 data values per sample per chip).
##
## Steven McKinney
## Molecular Oncology and Breast Cancer Program
## British Columbia Cancer Research Centre
## Oct 19, 2006
## smckinney at bccrc.ca
##
## Most credit goes to
## Natalie Thorne (snapCGH developer, currently in Dept of Oncology,
## University of Cambridge))
## Todd Richmond (Manager of Research Informatics, NimbleGen)
## Natalie Thorne and Todd Richmond provided most of the functions
## to allow us to get up and running in R for NimbleGen data.
## Many thanks to them.
##
## Note: NimbleGen's business model in part aims to stay flexible, to
## allow customization of chip features when necessary. NimbleGen
## therefore does not publish data or software specifications, and
## reserves the right to change file and data formats. It is thus
## difficult to produce code similar to that available for the
## affy chips, as NimbleGen specifications can change.
## You should thus carefully check that your data formats
## work with this code, and modify these scripts and functions
## accordingly.
## These scripts and functions come with no guarantee whatsoever.
## However, help is available - various Bioconductor and NimbleGen
## personnel certainly bent over backwards to help us out, for
## which we are most appreciative.
###---------------------------------------------------------------------
## Load required R/BioConductor libraries
library(snapCGH)
library(limma)
library(IDPmisc)
## Assuming that NimblegenUtilityFunctions.R is in working
## directory. Add appropriate file path otherwise.
source("NimblegenUtilityFunctions.R")
## NimRootDir should point to the directory containing
## all the NimbleGen data material as delivered on CD
## or DVD. This root directory will contain a sample key
## file, e.g. ORDER_ID3051_SampleKey.txt, and a subdirectory
## "Raw_data_files" containing the 'pair' data files with names
## such as "63392_532.pair".
## You may need to specify a different path here, if your
## NimbleGen data is not in the subdirectory "Raw_data_files"
## of the R working directory.
NimRootDir <- getwd()
## readSampleKey will by default look for a file whose name
## contains "SampleKey.txt", e.g. ORDER_ID3051_SampleKey.txt.
## If your NimbleGen order comes split over several CDs, you may have
## to manually create one SampleKey file containing all the Chip
## information for the chip data you want to process.
samplekey <- readSampleKey(path = NimRootDir, row.names = NULL)
## NimblegenTargetFiles builds the names of the data files
## based on the information from the SampleKey file.
data.files <-
NimblegenTargetFiles(samplekey,
path = paste(NimRootDir,
"/Raw_Data_Files",
sep=""))
## NimblegenTargets retrieves sample descriptive information, e.g.
## sample type and fluorescent labeling.
targets <- as.data.frame(NimblegenTargets(samplekey, path = NimRootDir))
## read.nimblegen() reads in all the 'pair' raw data intensity files.
RG <- read.nimblegen(files = data.files)
## Some of our order was customized, and the returned data files had
## extra quality control spot intensities in the raw data files.
## We want to examine only the spots that map to chromosomes.
## Data for these spots have PROBE_ID labels that start with "CHR",
## so we keep these and filter out the rest.
if ( length(chrIDxs <- grep("^CHR", as.character(RG$genes$PROBE_ID))) < nrow(RG) ) {
RG <- RG[chrIDxs, ]
}
if ( exists("chrIDxs") ) rm(chrIDxs)
## Our data is from female samples, so we filter out Y chromosome
## data as well.
if( length(chrYIDxs <- grep("^CHRY", as.character(RG$genes$PROBE_ID))) ) {
RG <- RG[-chrYIDxs, ]
}
if ( exists("chrYIDxs") ) rm(chrYIDxs)
## After filtering, we need to get rid of unnecessary factor classes
## for several of the variables.
RG$genes$SEQ_ID <- as.factor(as.character(RG$genes$SEQ_ID))
RG$genes$GENE_EXPR_OPTION <- as.factor(as.character(RG$genes$GENE_EXPR_OPTION))
RG$genes$PROBE_ID <- as.character(RG$genes$PROBE_ID)
## Produce heatmaps of the raw data.
NimblegenHeatmaps(RG)
## Produce spotmaps of the raw data to examine spot locations
## and chip features.
NimblegenHeatmaps(input = RG, filesuffix = "spotmap.pdf", ncolors = 1)
## Now we are ready to start processing the raw data.
## readPositionalInfo() sets up probe position information.
RG <- readPositionalInfo(RG, source = "nimblegen")
## Now we are ready to normalize. NimbleGen uses qspline, and we confirmed
## that no background subtraction plus within-chip q-spline normalization
## works well for our situation.
## Set up a matrix to hold normalization data.
pmi.normmtx <- matrix(NA, nrow = nrow(RG$G), ncol = 2)
for (iptr in seq(ncol(RG$G))) {
## Perform q-spline normalization for each chip.
pmi.normmtx[] <- normalize.qspline(as.matrix(cbind(RG$G[, iptr], RG$R[, iptr])))
RG$G[, iptr] <- pmi.normmtx[, 1]
RG$R[, iptr] <- pmi.normmtx[, 2]
}
rm(pmi.normmtx)
## Set up the 'design' component, as discussed in snapCGH
## documentation. In our case, CY5 (Red) is the reference sample.
RG$design <- rep(-1, nrow(targets))
## Get MA form for data
MA <- MA.RG(RG)
## MA plots
## plotMA(): limma library standard MA plot routine.
## Our chips have about 380000 data values, so standard
## MA plots are mostly black. Graphs and output files
## are large. Adobe reader chokes on these standard plotMA plots.
pdf(file = "plotMA.pdf")
lapply(seq(ncol(MA)),
function(x) { plotMA(MA, array = x)
mtext("MA after q-spline normalization") })
dev.off()
## I found MA 'heatmap' scatterplots to work well for this
## plot conumdrum. I've set up alternative MA plot routines
## using ixyplot() from package IDPmisc.
## iplot() produces an ixyplot with heat legend. I'd like
## to have vertical and horizontal lines in the plot, just haven't
## had the time to work them in.
pdf(file = "iplotMA.pdf")
lapply(seq(ncol(MA)), function(x) iplotMA(MA, array = x))
dev.off()
## To get e.g. horizontal line at 0, I use ImageplotMA for now.
pdf(file = "ImageplotMA.pdf")
lapply(seq(ncol(MA)), function(x) {ImageplotMA(MA, array = x); abline(h = 0)})
dev.off()
## Set up the data (e.g. put the clone data in order along the chromosomes).
MA.process <- processCGH(MA,
maxChromThreshold = 23, minChromThreshold = 1,
method.of.averaging = NULL, ID = "PROBE_ID")
## Run segmentation analysis.
SegInfo.HomHMM <- runHomHMM(MA.process, criteria = "AIC")
###---------------------------------------------------------------------
## End processNimblegenData.R
###---------------------------------------------------------------------
###---------------------------------------------------------------------
## NimblegenUtilityFunctions.R
## Utility functions for reading and processing NimbleGen data
## from NimbleGen pair files.
##
## Steven McKinney
## Molecular Oncology and Breast Cancer Program
## British Columbia Cancer Research Centre
## Oct 19, 2006
## smckinney at bccrc.ca
##
## Most credit goes to
## Natalie Thorne (snapCGH developer, currently in Dept of Oncology,
## University of Cambridge))
## Todd Richmond (Manager of Research Informatics, NimbleGen)
## Natalie Thorne and Todd Richmond provided most of the functions
## to allow us to get up and running in R for NimbleGen data.
## Many thanks to them.
##
## Note: NimbleGen's business model in part aims to stay flexible, to
## allow customization of chip features when necessary. NimbleGen
## therefore does not publish data or software specifications, and
## reserves the right to change file and data formats. It is thus
## difficult to produce code similar to that available for the
## affy chips, as NimbleGen specifications can change.
## You should thus carefully check that your data formats
## work with this code, and modify these scripts and functions
## accordingly.
## These scripts and functions come with no guarantee whatsoever.
## However, help is available - various Bioconductor and NimbleGen
## personnel certainly bent over backwards to help us out, for
## which we are most appreciative.
###---------------------------------------------------------------------
"readSampleKey" <-
function(file = "SampleKey.txt",
path = NULL,
sep = "\t",
row.names = "SAMPLE_LABEL",
quote = "\"",
...)
{
## Find the sample key file.
## Nimblegen file naming convention is of the form
## ORDER_ID3051_SampleKey.txt
if (is.null(path)) path <- getwd()
## Get a list of files in the <path> directory.
filenames.vec <- list.files(path)
## Use <file> as a mask
SampleKeyPosn <- grep(file, filenames.vec)
if (length(SampleKeyPosn) == 0) {
msg <- paste("\nNo sample key file found containing string '",
file,
"'",
"\nin directory\n",
path,
"\n\n",
sep = "")
stop(msg)
}
if (length(SampleKeyPosn) > 1) {
msg <- paste("\nMultiple sample key files found containing string '",
file,
"'",
"\nin directory\n",
path,
"\n\n",
sep = "")
stop(msg)
}
file <- filenames.vec[SampleKeyPosn]
if (!is.null(path)) file <- file.path(path, file)
SampleKeyHeader <- readSampleKeyHeader(file)
if(is.null(row.names)) {
## Need to handle different experiment scenario where
## same reference sample is used in copy number comparisons.
tab <- read.table(file, header = TRUE, as.is = TRUE, sep = sep,
quote = quote, fill = TRUE, ...)
} else {
## Natalie Thorne's original read table for her style data.
row.num <- which(row.names == SampleKeyHeader)
tab <- read.table(file, header = TRUE, as.is = TRUE, sep = sep,
quote = quote, fill = TRUE, row.names = row.num,...)
##row.names(tab) <- removeExt(tab[,row.num])
}
tab
}
"readSampleKeyHeader" <- function(file)
{
con <- file(file, "r")
on.exit(close(con))
out <- list()
i <- 0
txt <- readLines(con, n = 1)
if (txt == "")
stop("Not SampleKey data file: input stopped at blank line")
as.vector(unlist(strsplit(txt,split="\t")))
}
"NimblegenTargetFiles" <- function(samplekey, path = NULL) {
chipid <- unique(samplekey[,which("CHIP_ID" == colnames(samplekey))])
cy3 <- paste(chipid,"532",sep="_")
cy3 <- paste(cy3,"pair",sep=".")
cy5 <- paste(chipid,"635",sep="_")
cy5 <- paste(cy5,"pair",sep=".")
if(!is.null(path)){
cy3 <- paste(path,cy3,sep="/")
cy5 <- paste(path,cy5,sep="/")
}
cbind(cy3,cy5)
}
"NimblegenTargets" <-
function(samplekey,
path = NULL)
{
## Get sample information from sample key file.
## CHIP_ID: All should occur in pairs.
unique.chipids <- unique(samplekey[, "CHIP_ID"])
if ( !all(table(samplekey[, "CHIP_ID"]) == 2) ) {
msg <- paste("The following chip IDs do not occur in pairs:\n",
labels(which(table(samplekey[, "CHIP_ID"]) != 2)),
"\n",
sep = "")
stop(msg)
}
n.chips <- length(unique.chipids)
out <- matrix("", ncol = 5, nrow = n.chips)
colnames(out) <- c(
"Array",
"OID",
"Tissue",
"Cy3",
"Cy5"
)
for ( i in seq(along = unique.chipids)) {
chipid.i <- unique.chipids[i]
chipid.i.rows <- which(samplekey[, "CHIP_ID"] == chipid.i)
out[i, "Array"] <- chipid.i
out[i, "OID"] <- samplekey[chipid.i.rows[1], "ORD_ID"]
out[i, "Tissue"] <-
unlist(strsplit(samplekey[chipid.i.rows[1], "SAMPLE_DESCRIPTION"], "-"))[1]
out[i, "Cy3"] <-
samplekey[which(samplekey[, "CHIP_ID"] == chipid.i &
tolower(samplekey[, "DYE"]) == "cy3"),
"SAMPLE_DESCRIPTION"]
out[i, "Cy5"] <-
samplekey[which(samplekey[, "CHIP_ID"] == chipid.i &
tolower(samplekey[, "DYE"]) == "cy5"),
"SAMPLE_DESCRIPTION"]
}
return(out)
}
"read.nimblegen" <-
function (files,
path = NULL,
ext = NULL,
names = NULL,
columns = NULL,
wt.fun = NULL,
verbose = TRUE,
sep = "\t",
quote = "\"",
...)
{
if (is.null(dim(files))) {
if (length(files)%%2 == 0)
files <- matrix(files, ncol = 2, byrow = TRUE)
else stop("Odd number of files: should be two data files for each array")
}
else {
files <- as.matrix(files)
if (ncol(files) != 2)
stop("Need a two column matrix of file names")
}
if (!is.null(ext))
files <- array(paste(files, ext, sep = "."), dim(files))
narrays <- nrow(files)
if (is.null(columns))
columns <- list(f = "PM")
if (is.null(names))
names <- removeExt(basename(files[, 1]))
fullname <- files[1, 1]
if (!is.null(path))
fullname <- file.path(path, fullname)
skip <- 1
obj <- read.table(fullname, skip = skip, header = TRUE,
sep = sep, quote = quote, check.names = FALSE, comment.char = "",
fill = TRUE,...)
nspots <- dim(obj)[1]
Y <- matrix(0, nspots, narrays)
colnames(Y) <- names
RG <- list(R = Y, G = Y, Rb = Y, Gb = Y)
if (!is.null(wt.fun))
RG$weights <- Y
for (i in 1:narrays) {
fullname <- files[i, 1]
if (!is.null(path))
fullname <- file.path(path, fullname)
obj <- read.table(fullname, skip = skip, header = TRUE,
sep = sep, quote = quote, check.names = FALSE, comment.char = "",
fill = TRUE,...)
if (verbose)
cat(paste("Read", fullname, "\n"))
if (i == 1)
RG$genes <- obj[, c("IMAGE_ID", "GENE_EXPR_OPTION", "SEQ_ID",
"PROBE_ID", "POSITION", "X","Y","MATCH_INDEX","SEQ_URL")]
## Ensure that MATCH_INDEX alignment occurs for all files.
if(!all.equal(RG$genes$MATCH_INDEX, obj[, "MATCH_INDEX"]))
stop(paste("MATCH_INDEX out of sequence for file\n", fullname))
RG$G[, i] <- obj[, columns$f]
fullname <- files[i, 2]
if (!is.null(path))
fullname <- file.path(path, fullname)
obj <- read.table(fullname, skip = skip, header = TRUE,
sep = sep, quote = quote, check.names = FALSE, comment.char = "",
fill = TRUE, ...)
if (verbose)
cat(paste("Read", fullname, "\n"))
## Ensure that MATCH_INDEX alignment occurs for all files.
if(!all.equal(RG$genes$MATCH_INDEX, obj[, "MATCH_INDEX"]))
stop(paste("MATCH_INDEX out of sequence for file\n", fullname))
RG$R[, i] <- obj[, columns$f]
if (!is.null(wt.fun))
RG$weights[, i] <- wt.fun(obj)
}
new("RGList", RG)
}
"NimblegenHeatmaps" <-
function(input,
filesuffix = "heatmap.pdf",
ngAspectRatio = 1.024,
ncolors = 123
)
{
## NimblegenHeatmaps: Produce heat map for each R and G component
## of input.
## Usage: NimblegenHeatmaps(input = RG, filesuffix = "heatmap.pdf", ngAspectRatio = 1.024)
## NimblegenHeatmaps(RG)
## Args:
## input: RGList data object
## filesuffix: file name suffix for output pdf files.
## ngAspectRatio: aspect ratio to produce plots with individual spots
## appearing as squares.
## Verify input:
if ( class(input) != "RGList") stop("input must be an RGList object.")
## Re-orient the data to match NimbleGen orientation
idxmat <- cbind(input$genes$X, max(input$genes$Y) + 1 - input$genes$Y)
## Set up zmat with a border of 2 empty rows/cols
## so heatmap does not print on border box.
zmat <- matrix(as.numeric(NA),
nrow = max(input$genes$X) + 4,
ncol = max(input$genes$Y) + 4)
## Raw data heat maps
## in NimbleGen orientation
on.exit(cat("\n"))
cat("\nProcessing heat maps...\n")
for (iptr in seq(ncol(input$G))) {
chipID <- dimnames(input$G)[[2]][iptr]
filename <- paste("ChipID_", chipID, "_G_", filesuffix, sep = "")
pdf(file = filename)
cat("G:ChipNo = ", iptr, "(ChipID = ", chipID, ") ")
zmat[idxmat + 2] <- input$G[, iptr]
## colors
if(ncolors == 1) {
low <- col2rgb("green")/255
} else {
low <- col2rgb("white")/255
}
high <- col2rgb("green")/255
col <- rgb(seq(low[1], high[1], len = ncolors),
seq(low[2], high[2], len = ncolors),
seq(low[3], high[3], len = ncolors))
image(x = (-1:770), y = (-1:1026), z = zmat,
col = col, xlab = "X", ylab = "Y", asp = ngAspectRatio)
title(main = paste("G Sample = ", chipID))
dev.off()
}
cat("\n")
for (iptr in seq(ncol(input$R))) {
chipID <- dimnames(input$R)[[2]][iptr]
filename <- paste("ChipID_", chipID, "_R_", filesuffix, sep = "")
pdf(file = filename)
cat("R:ChipNo = ", iptr, "(ChipID = ", chipID, ") ")
zmat[idxmat + 2] <- input$R[, iptr]
## colors
if(ncolors == 1) {
low <- col2rgb("red")/255
} else {
low <- col2rgb("white")/255
}
high <- col2rgb("red")/255
col <- rgb(seq(low[1], high[1], len = ncolors),
seq(low[2], high[2], len = ncolors),
seq(low[3], high[3], len = ncolors))
image(x = (-1:770), y = (-1:1026), z = zmat,
col = col, xlab = "X", ylab = "Y", asp = ngAspectRatio)
title(main = paste("R Sample = ", chipID))
dev.off()
}
}
###----------------------------------------------
iplotMA <-
function (MA,
array = 1,
xlab = "A",
ylab = "M",
main = colnames(MA)[array],
xlim = NULL,
ylim = NULL,
status,
values,
pch,
col,
cex,
legend = TRUE,
zero.weights = FALSE,
...)
{
## Plot function for MA plots with many rows.
## limma library plotMA function yields largely black
## graphs that are difficult to render in pdf readers,
## printers and documents when hundreds of thousands
## of data values are plotted.
## Need ixyplot() from package IDPmisc.
require(IDPmisc)
switch(class(MA),
RGList = {
MA <- MA.RG(MA[, array])
array <- 1
x <- MA$A
y <- MA$M
w <- MA$w
},
MAList = {
x <- as.matrix(MA$A)[, array]
y <- as.matrix(MA$M)[, array]
if (is.null(MA$weights))
w <- NULL
else w <- as.matrix(MA$weights)[, array]
},
list = {
if (is.null(MA$A) || is.null(MA$M))
stop("No data to plot")
x <- as.matrix(MA$A)[, array]
y <- as.matrix(MA$M)[, array]
if (is.null(MA$weights))
w <- NULL
else w <- as.matrix(MA$weights)[, array]
},
MArrayLM = {
x <- MA$Amean
y <- as.matrix(MA$coefficients)[, array]
if (is.null(MA$weights))
w <- NULL
else w <- as.matrix(MA$weights)[, array]
},
matrix = {
narrays <- ncol(MA)
if (narrays < 2)
stop("Need at least two arrays")
if (narrays > 5)
x <- apply(MA, 1, median, na.rm = TRUE)
else x <- rowMeans(MA, na.rm = TRUE)
y <- MA[, array] - x
w <- NULL
},
exprSet = {
narrays <- ncol(MA at exprs)
if (narrays < 2)
stop("Need at least two arrays")
if (narrays > 5)
x <- apply(MA at exprs, 1, median, na.rm = TRUE)
else x <- rowMeans(MA at exprs, na.rm = TRUE)
y <- MA at exprs[, array] - x
w <- NULL
if (missing(main))
main <- colnames(MA at exprs)[array]
},
stop("MA is invalid object")
)
if (!is.null(w) && !zero.weights) {
i <- is.na(w) | (w <= 0)
y[i] <- NA
}
if (is.null(xlim))
xlim <- range(x, na.rm = TRUE)
if (is.null(ylim))
ylim <- range(y, na.rm = TRUE)
if (missing(status))
status <- MA$genes$Status
ixyplot(x, y, xlab = xlab, ylab = ylab, main = main)
## abline(h = 0) ## Will have to put this in ixyplot before key legend plot
invisible()
}
ImageplotMA <-
function (MA, array = 1, xlab = "A", ylab = "M", main = colnames(MA)[array],
xlim = NULL, ylim = NULL, status, values, pch, col, cex,
legend = TRUE, zero.weights = FALSE, ...)
{
## Plot function for MA plots with many rows.
## limma library plotMA function yields largely black
## graphs that are difficult to render in pdf readers,
## printers and documents when hundreds of thousands
## of data values are plotted.
## Need ixyplot() from package IDPmisc.
require(IDPmisc)
switch(class(MA), RGList = {
MA <- MA.RG(MA[, array])
array <- 1
x <- MA$A
y <- MA$M
w <- MA$w
}, MAList = {
x <- as.matrix(MA$A)[, array]
y <- as.matrix(MA$M)[, array]
if (is.null(MA$weights))
w <- NULL
else w <- as.matrix(MA$weights)[, array]
}, list = {
if (is.null(MA$A) || is.null(MA$M))
stop("No data to plot")
x <- as.matrix(MA$A)[, array]
y <- as.matrix(MA$M)[, array]
if (is.null(MA$weights))
w <- NULL
else w <- as.matrix(MA$weights)[, array]
}, MArrayLM = {
x <- MA$Amean
y <- as.matrix(MA$coefficients)[, array]
if (is.null(MA$weights))
w <- NULL
else w <- as.matrix(MA$weights)[, array]
}, matrix = {
narrays <- ncol(MA)
if (narrays < 2)
stop("Need at least two arrays")
if (narrays > 5)
x <- apply(MA, 1, median, na.rm = TRUE)
else x <- rowMeans(MA, na.rm = TRUE)
y <- MA[, array] - x
w <- NULL
}, exprSet = {
narrays <- ncol(MA at exprs)
if (narrays < 2)
stop("Need at least two arrays")
if (narrays > 5)
x <- apply(MA at exprs, 1, median, na.rm = TRUE)
else x <- rowMeans(MA at exprs, na.rm = TRUE)
y <- MA at exprs[, array] - x
w <- NULL
if (missing(main))
main <- colnames(MA at exprs)[array]
}, stop("MA is invalid object"))
if (!is.null(w) && !zero.weights) {
i <- is.na(w) | (w <= 0)
y[i] <- NA
}
if (is.null(xlim))
xlim <- range(x, na.rm = TRUE)
if (is.null(ylim))
ylim <- range(y, na.rm = TRUE)
if (missing(status))
status <- MA$genes$Status
plot(x, y, xlab = xlab, ylab = ylab, main = main, xlim = xlim,
ylim = ylim, type = "n", ...)
if (is.null(status) || all(is.na(status))) {
if (missing(pch))
pch = 16
if (missing(cex))
cex = 0.3
Image(x, y)
}
else {
if (missing(values)) {
if (is.null(attr(status, "values")))
values <- names(sort(table(status), decreasing = TRUE))
else values <- attr(status, "values")
}
sel <- !(status %in% values)
nonhi <- any(sel)
if (nonhi)
Image(x[sel], y[sel])
nvalues <- length(values)
if (missing(pch)) {
if (is.null(attr(status, "pch")))
pch <- rep(16, nvalues)
else pch <- attr(status, "pch")
}
if (missing(cex)) {
if (is.null(attr(status, "cex"))) {
cex <- rep(1, nvalues)
if (!nonhi)
cex[1] <- 0.3
}
else cex <- attr(status, "cex")
}
if (missing(col)) {
if (is.null(attr(status, "col"))) {
col <- nonhi + 1:nvalues
}
else col <- attr(status, "col")
}
pch <- rep(pch, length = nvalues)
col <- rep(col, length = nvalues)
cex <- rep(cex, length = nvalues)
for (i in 1:nvalues) {
sel <- status == values[i]
Image(x[sel], y[sel])
}
if (legend) {
if (is.list(pch))
legend(x = xlim[1], y = ylim[2], legend = values,
fill = col, col = col, cex = 0.9)
else legend(x = xlim[1], y = ylim[2], legend = values,
pch = pch, , col = col, cex = 0.9)
}
}
invisible()
}
######################################################
###---------------------------------------------------------------------
## End NimblegenUtilityFunctions.R
###---------------------------------------------------------------------
-------------- next part --------------
An embedded and charset-unspecified text was scrubbed...
Name: processNimblegenData.R
Url: https://stat.ethz.ch/pipermail/bioconductor/attachments/20061019/669af712/attachment.pl
-------------- next part --------------
An embedded and charset-unspecified text was scrubbed...
Name: NimblegenUtilityFunctions.R
Url: https://stat.ethz.ch/pipermail/bioconductor/attachments/20061019/669af712/attachment-0001.pl
More information about the Bioconductor
mailing list