Hi Nathalie -- On 06/03/2011 08:52 AM, Nathalie Conte wrote: > HI, I a new in R and bioconductor and need somebody' s input on these > issues ( see below, 2 questions), I am sorry about the rather long It's always daunting to jump in to the middle of someone's problem; it helps to have a working, real example. I took a look at the example in ?segment and came up with library(DNAcopy) set.seed(25) genomdat <- rnorm(500, sd=0.1) + rep(c(-0.2,0.1,1,-0.5,0.2,-0.5,0.1,-0.2), c(137,87,17,49,29,52,87,42)) chrom <- rep(1:2, c(290,210)) maploc <- c(1:290, 1:210) plot(segment(CNA(genomdat, chrom, maploc))) Is that the right starting point? some more below... > email, I am trying to put as much info as I can > first my system: > > sessionInfo() > R version 2.11.1 (2010-05-31) > x86_64-unknown-linux-gnu > > locale: > [1] LC_CTYPE=en_GB.UTF-8 LC_NUMERIC=C > [3] LC_TIME=en_GB.UTF-8 LC_COLLATE=en_GB.UTF-8 > [5] LC_MONETARY=C LC_MESSAGES=C > [7] LC_PAPER=en_GB.UTF-8 LC_NAME=C > [9] LC_ADDRESS=C LC_TELEPHONE=C > [11] LC_MEASUREMENT=en_GB.UTF-8 LC_IDENTIFICATION=C > > attached base packages: > [1] stats graphics grDevices utils datasets methods base > > other attached packages: > [1] DNAcopy_1.24.0 > > loaded via a namespace (and not attached): > [1] tools_2.11.1 > > > I am trying to plot some cgharray data (mouse agilent 244) and I need to > make plots. > I went for the dNAcopy package, following normalisation using limma, I > am first creating a copy number array using CNA function from DNAcopy, I > have used this script (based on the cghMCR package pdf > instructions).Ma_betAr is my normalised data. > >chrom <- gsub("chr([0-9XY]+):.*", "\\1", Ma_BetAr$genes[, > "SystematicName"]) probably 'sub' instead of 'gsub' > >dropMe <- c(which(!chrom %in% c(1:19, "X", "Y")), > which(Ma_BetAr$genes[, "ControlType"] != 0)) > > >cna <- CNA(Ma_BetAr$M[-dropMe, ], > gsub("chr([0-9XY]+):.*", "\\1", Ma_BetAr$genes[-dropMe, "SystematicName"]), > as.numeric(gsub(".*:([0-9]+)-.*", "\\1", > Ma_BetAr$genes[-dropMe, "SystematicName"])), > data.type = "logratio", sampleid = colnames(Ma_BetAr$M)) > > >smooth.cna<-smooth.CNA(cna) smoothing the cna > then this smoothed CNA is segment using the segment function > >segData <- segment(smooth.CNA(cna)) > > then this segData is plotted using plot(segData, plot.type = "w") . > At this point there is a plot(log2R in Yaxis and chromosome in X) which > is created that will alternate the chromosomes alternate in different > colours and are sorted in this order: 1 then 10 then 11, 12, 13 ,14 ,15, > 16, 17 ,18 ,19 ,2 ,3 ,4,... > Question1:I would rather have them in the natural order: 1,2,3,4.....Is > there an easy way to change the plot order? I saw at the CNA step on the > manual(DNAcopy) there is a mention of natural order but I don't quite > understand what needs to be done:see below > Usage > CNA(genomdat, chrom, maploc, data.type=c("logratio","binary"), > sampleid=NULL) > ## S3 method for class 'CNA': > print(x, ...) > Arguments > a vector or matrix of data from array-CGH, ROMA, or other copy number ex- > genomdat > periments. If it is a matrix the rows correspond to the markers and the > columns > to the samples. > the chromosomes (or other group identifier) from which the markers came. > Vec- > chrom > tor of length same as the number of rows of genomdat. If one wants the > chro- > mosomes to be ordered in the natural order, this variable should be > numeric or > ordered category. In my example, I see that chrom <- rep(c("chr11", "chr2"), c(290,210)) plot(segment(CNA(genomdat, chrom, maploc)), plot.type="whole") gives me the wrong order (because "chr11" comes before "chr2" if I 'sort' these). I suspect the author of the package intended for something different, but if I create a map between my chromosome names and the order in which I'd like them to appear... map <- 1:24 names(map) <- paste("chr", c(1:22, "X", "Y"), sep="") and then recode my chromosome names to their integer representation chrom <- map[rep(c("chr11", "chr2"), c(290,210))] plot(segment(CNA(genomdat, chrom, maploc)), plot.type="whole") everyone is happy. I think the author intended that it would possible to create a (ordered?) factor, where the levels indicate the sort order. chrom <- factor(rep(c("chr11", "chr2"), c(290,210)), levels=paste("chr", c(1:22, "X", "Y"), sep="")) but somewhere things are going awry. > the locations of marker on the genome. Vector of length same as the > number of > maploc > rows of genomdat. This has to be numeric. > logratio (aCGH, ROMA, etc.) or binary (LOH). > data.type > sample identifier. If missing the samples are named by prefixing > "Sample" to > sampleid > consecutive integers. > object returned by CNA > x > arguments to be passed onto print command called within. > ... > > > Question2: Is it possible at this point to plot the probe using a colour > that will correspond to their log2R, example if a probe log2R is more > than 0.2 plot it in red using a script? I don't think so using plot.DNAcopy. But in some ways a little creativity goes a long way (this is probably too much for an answer, but maybe inspires...) library(RColorBrewer) # nice palettes; see http://colorbrewer2.org/ library(lattice) pal <- brewer.pal(9, "OrRd")[-1] # drop lightest color breaks <- with(seg$data, { seq(min(Sample.1), max(Sample.1), length.out=length(pal) + 1) }) xyplot(Sample.1 ~ maploc | as.factor(chrom), seg$data, panel=function(x, y, ..., subscripts, segments, pal, breaks) { col <- pal[cut(y, breaks=breaks, labels=FALSE)] panel.xyplot(x, y, ..., col=col) currChrom <- unique(seg$data$chrom[subscripts]) segs <- segments[as.factor(segments$chrom) == currChrom,, drop=FALSE] for (i in seq_len(nrow(segs))) with(segs[i,,drop=FALSE], { llines(c(loc.start, loc.end), seg.mean, col="black") }) }, segments=seg$output, pch=20, pal=pal, breaks=breaks, layout=c(1, 2), strip=FALSE, strip.left=TRUE) Martin > > Many thanks in advance for any help > Nathalie > > > > > > -- Computational Biology Fred Hutchinson Cancer Research Center 1100 Fairview Ave. N. PO Box 19024 Seattle, WA 98109 Location: M1-B861 Telephone: 206 667-2793