Hi Henrik, Thanks a lot for the help! Sorry I have more questions. I am following "How to: Calculate total copy number ratios from total (non-polymorphic) signals" and "Vignette: Total copy-number segmentation (non-paired CBS)", but I am not sure if I do it correctly. I have two SNP6 datasets Tumor and HapMap270 and I want to use HpaMap270 as reference to go all the way to CBS step. So I do the following steps respectively. > ds1 <- doCRMAv2("HapMap270", chipType="GenomeWideSNP_6,Full") > ds2 <- doCRMAv2("Tumor", chipType="GenomeWideSNP_6,Full") After that, I do > dataSet <- "HapMap270" > tags <- "ACC,ra,-XY,BPN,-XY,AVG,A+B,FLN,-XY" > chipType <- "GenomeWideSNP_6" > dsN <- AromaUnitTotalCnBinarySet$byName(dataSet, tags=tags, chipType=chipType) > dataSet <- "Tumor" > tags <- "ACC,ra,-XY,BPN,-XY,AVG,A+B,FLN,-XY" > chipType <- "GenomeWideSNP_6" > dsT <- AromaUnitTotalCnBinarySet$byName(dataSet, tags=tags, chipType=chipType) > dfR <- getAverageFile(dsN) # ? > dsTR <- exportTotalCnRatioSet(dsT, ref=dfR) # ? Would the above two steps work? My question is how to go ahead from here to do CBS and will sm <- CbsModel(dsTR) work? Thanks again for your help! Sean
On Sun, Jul 28, 2013 at 4:28 PM, Henrik Bengtsson < henrik.bengts...@aroma-project.org> wrote: > Hi. > > On Fri, Jul 26, 2013 at 8:02 AM, sean nj <njs...@gmail.com> wrote: > > Hi guys, > > > > I have a question regarding how to calculate raw copy numbers using > common > > reference instead of average of all samples of the study. Basically I > want > > to use average of HapMap270 samples as reference for all further copy > number > > calculations. > > > > I have a bunch HapMap270 snp6 cel files and I followed Vignette: > Estimation > > of total copy numbers using the CRMA v2 method (10K-CytoScanHD) to Step > 5 - > > Calculation of raw copy numbers, and generated ceR and saved it as a > RData > > file ceR.Rdata. > > It's important to understand that almost all objects in the Aroma > framework are basically "pointers" to external files. For instance, > your 'ceR', which I assume you've got from something like ceR <- > getAverageFile(ces), is referring to the file with pathname > getPathname(ceR). More below... > > > > > My first question is, how to use this data for any future copy number > > analysis? My guess is that instead of calculating the ceR from the sample > > set I can just load the ceR.RData file I saved and use it. Right? > > First of all, please note that when do ceR <- getAverageFile(ces) on > the same data set 'ces', the result is already available on file and > it will be quickly found and returned. In other words, it will not > recalcuate the averages again [unless you do ceR <- > getAverageFile(ces, force=TRUE)]. > > However, I do understand that you may not want to have to keep a large > 'ces' data set around, when you're only interested in the pooled > average. In that case, I would copy the file containing the "average" > to a new data set. Currently, this is not straightforward in Aroma > (I'll think about something), but you can do the following: > > # Calculate the pooled average > > ceR <- getAverageFile(cesN); > > # Copy this file to plmData/HapMap270,pooled/GenomeWideSNP_6/, e.g. > > filename <- getFilename(ceR); > > filename > [1] ".average-intensities-median-mad,d03faaf8b707a97c4e43381b1a5d1ef2.CEL" > > rootPath <- getParent(getPath(cesN), depth=2L); > > dataSet <- "HapMap270,pooled"; > > chipType <- getChipType(ceR, fullname=FALSE); > > path <- file.path(rootPath, dataSet, chipType); > > path > [1] "plmData/HapMap270,pooled/GenomeWideSNP_6" > > mkdirs(path); > > copyFile(getPathname(ceR), file.path(path, filename)); > > With this done, you can then grab this pooled reference as: > > > library("aroma.affymetrix") > > path <- "plmData/HapMap270,pooled/GenomeWideSNP_6"; > > filename <- > ".average-intensities-median-mad,d03faaf8b707a97c4e43381b1a5d1ef2.CEL"; > > ceR <- CnChipEffectFile(filename, path=path); > > Note, when you save 'ceR', you are basically saving the reference to > the file. Yes, you can load it later, but make sure not to move it, > otherwise you'll get some type of "file not found" error. > > > saveObject(ceR, "HapMap270,GenomeWideSNP_6,reference.Rdata"); > > If already saved, and file not moved, you can then do: > > > library("aroma.affymetrix"); > > ceR <- loadObject("HapMap270,GenomeWideSNP_6,reference.Rdata"); > > All this is very ad hoc (=non-aroma style), and as I said, I'll see if > I can come up with a cleaner solution for storing and retrieving > pooled averages. > > > > > My second question is, how to go ahead from there to calculate the > relative > > copy numbers for all unit from all samples? The two examples given in the > > Vignette are for one unit from one sample and for a few unit on > chromosome 2 > > for one sample. What is the function to retrieve all units on all > > chromosomes instead of units <- getUnitsOnChromosome(gi, chromosome=2, > > region=c(81,86)*1e6)? > > You can set 'units' to NULL to retrieve all loci, i.e. no need to use > getUnitsOnChromosome(). FYI, units <- NULL will give the same data as > with units <- 1:nbrOfUnits(gi). > > > And what is the function to retrieve all samples > > instead of ce <- getFile(cesN, indexOf(cesN, "NA06985"))? > > Hmm... not clear what you mean. All samples are in 'cesN', and you do > need to iterate over them somehow. Is this what you're looking for? > > for (ii in seq_along(cesN)) { > ce <- getFile(cesN, ii) > ... > } > > Or are you asking how to extract the data from all samples? Then you can > do: > > theta <- extractTheta(cesN, units=units) > > but be careful because that loads a lot of data into memory. > > Hope this helps, > > Henrik > > > > > Thanks a lot for the help, > > > > Sean > > > > -- > > -- > > When reporting problems on aroma.affymetrix, make sure 1) to run the > latest > > version of the package, 2) to report the output of sessionInfo() and > > traceback(), and 3) to post a complete code example. > > > > > > You received this message because you are subscribed to the Google Groups > > "aroma.affymetrix" group with website http://www.aroma-project.org/. > > To post to this group, send email to aroma-affymetrix@googlegroups.com > > To unsubscribe and other options, go to > http://www.aroma-project.org/forum/ > > > > --- > > You received this message because you are subscribed to the Google Groups > > "aroma.affymetrix" group. > > To unsubscribe from this group and stop receiving emails from it, send an > > email to aroma-affymetrix+unsubscr...@googlegroups.com. > > For more options, visit https://groups.google.com/groups/opt_out. > > > > > > -- > -- > When reporting problems on aroma.affymetrix, make sure 1) to run the > latest version of the package, 2) to report the output of sessionInfo() and > traceback(), and 3) to post a complete code example. > > > You received this message because you are subscribed to the Google Groups > "aroma.affymetrix" group with website http://www.aroma-project.org/. > To post to this group, send email to aroma-affymetrix@googlegroups.com > To unsubscribe and other options, go to > http://www.aroma-project.org/forum/ > > --- > You received this message because you are subscribed to the Google Groups > "aroma.affymetrix" group. > To unsubscribe from this group and stop receiving emails from it, send an > email to aroma-affymetrix+unsubscr...@googlegroups.com. > For more options, visit https://groups.google.com/groups/opt_out. > > > -- -- When reporting problems on aroma.affymetrix, make sure 1) to run the latest version of the package, 2) to report the output of sessionInfo() and traceback(), and 3) to post a complete code example. You received this message because you are subscribed to the Google Groups "aroma.affymetrix" group with website http://www.aroma-project.org/. To post to this group, send email to aroma-affymetrix@googlegroups.com To unsubscribe and other options, go to http://www.aroma-project.org/forum/ --- You received this message because you are subscribed to the Google Groups "aroma.affymetrix" group. To unsubscribe from this group and stop receiving emails from it, send an email to aroma-affymetrix+unsubscr...@googlegroups.com. For more options, visit https://groups.google.com/groups/opt_out.