Thanks a lot, Henrik! Cheers, Oscar Oscar M. Rueda, PhD. Postdoctoral Research Fellow, Caldas Lab, Breast Cancer Functional Genomics. University of Cambridge. Cancer Research UK Cambridge Institute. Li Ka Shing Centre, Robinson Way. Cambridge CB2 0RE England
Please note that my email address will be changing to: oscar.ru...@cruk.cam.ac.uk On 24/04/2013 10:25, "aroma-affymetrix@googlegroups.com" <aroma-affymetrix@googlegroups.com> wrote: > > Today's Topic Summary >Group: >http://groups.google.com/group/aroma-affymetrix/topics ><http://groups.google.com/group/aroma-affymetrix/topics> > >* PSCBS and callNTCN() <#group_thread_0> [1 Update] > >* callAB gives Error in bootstrapTCNandDHByRegion.PairedPSCBS ><#group_thread_1> [1 Update] > > > PSCBS and callNTCN() ><http://groups.google.com/group/aroma-affymetrix/t/94d01cba9523778c> > >Henrik Bengtsson <h...@biostat.ucsf.edu> > Apr 23 01:13PM -0700 > > > > >Oscar, sorry for the delay on this one. > > > > > >Your problem was somewhat related to Ingrid's non-paired data (see > > >other thread of today). I've verified that PSCBS v0.35.0 runs all the > > >way through with ROH, AB, LOH and NTCN calls using your data > > >('faultySample.Rbin'; 38,474,760 bytes). Update by: > > > > > >source(" >http://aroma-project.org/hbLite.R >"); > > >hbInstall("PSCBS"); > > > > > >For what's updated, see NEWS under PSCBS after doing help.start(). > > > > > >The script I used to verify your paired tumor-normal data is: > > > > > >stopifnot(packageVersion("PSCBS") >= "0.35.0"); > > >library("PSCBS"); > > >library("R.devices"); > > >devOptions("png", width=1024); > > >setOption("devEval/args/force", FALSE); > > > > > >sampleName <- "faultySample"; > > > > > >## Load the data you send me several weeks ago > > >pathname <- sprintf("%s.Rbin", sampleName); > > >fit <- loadObject(pathname); > > >data <- getLocusData(fit); > > >rm(fit); > > > > > >## Locate centromeres and other large gaps in data > > >gaps <- findLargeGaps(data, minLength=1e6); > > >knownSegments <- gapsToSegments(gaps); > > >## => dim(knownSegments) == c(64, 4) > > > > > >## Paired PSCBS segmentation > > >## 'betaT' is already normalized using TumorBoost => tbn=FALSE. > > >fit <- segmentByPairedPSCBS(data, knownSegments=knownSegments, > > >tbn=FALSE, avgDH="median", seed=4874, verbose=-10); > > >printf("Number of segments: %d\n", nbrOfSegments(fit)); > > >## Number of segments: 481 > > > > > >## Identify run of homozygosity (ROH) in germline > > >fit <- callROH(fit, verbose=-10); > > >printf("Number of ROH segments: %d\n", sum(getSegments(fit)$rohCall, > > >na.rm=TRUE)); > > >## Number of ROH segments: 44 > > > > > >## Estimate global background level in [0,1] (due to normal > > >contamination and more) > > >kappa <- estimateKappa(fit, verbose=-10); > > >printf("Kappa: %g\n", kappa); > > >## Kappa: 0.700498 => Large background signal > > > > > >## Call allelic balance > > >## (a) Estimate DH threshold for calling AB > > >deltaAB <- estimateDeltaAB(fit, kappa=kappa); # If skipped, will be > > >done internally > > >printf("Delta_AB: %g\n", deltaAB); > > >## Delta_AB: 0.176208 > > > > > >## (b) Call AB based on bootstrapped segment DH levels > > >fit <- callAB(fit, delta=deltaAB, verbose=-10); > > >printf("Number of AB segments: %d\n", sum(getSegments(fit)$abCall, >na.rm=TRUE)); > > >## Number of AB segments: 321 > > > > > >## Call loss of heterozygosity (LOH) > > >## (a) Estimate C_1 threshold for calling LOH > > >deltaLOH <- estimateDeltaLOH(fit); # If skipped, will be done internally > > >printf("Delta_LOH: %g\n", deltaLOH); > > >## Delta_LOH: 0.75476 > > > > > >## (b) Call LOH based on bootstrapped segment C_1 levels > > >fit <- callLOH(fit, delta=deltaLOH, verbose=-10); > > >printf("Number of LOH segments: %d\n", sum(getSegments(fit)$lohCall, > > >na.rm=TRUE)); > > >## Number of LOH segments: 71 > > > > > >## Call NTCN > > >## (a) Estimate the threshold for calling neutral TCN segments > > >## By shrinking 'scale', more segments will be non-NTCN. > > >deltaCN <- estimateDeltaCN(fit, scale=1, kappa=kappa); > > >printf("Delta_CN: %g\n", deltaCN); > > >## Delta_CN: 0.149751 > > > > > >## (b) Call NTCN based on bootstrapped segment TCN levels > > >fit <- callNTCN(fit, delta=deltaCN, verbose=-20); > > >printf("Number of NTCN segments: %d\n", sum(getSegments(fit)$ntcnCall, > > >na.rm=TRUE)); > > >## Number of NTCN segments: 296 > > > > > >toPNG(sampleName, tags="tracks,avgDH=median,ROH+AB+LOH+NTCN", >aspectRatio=0.6, { > > >plotTracks(fit); > > >}) > > > > > >/Henrik > > > > > > > > > > > > callAB gives Error in bootstrapTCNandDHByRegion.PairedPSCBS ><http://groups.google.com/group/aroma-affymetrix/t/e0d930961ad43ad2> > >Henrik Bengtsson <henrik.bengts...@aroma-project.org> > Apr 23 12:57PM -0700 > > > > >Hi, > > > > > >I could reproduce the error with the data you've sent me offline. > > >I've located the bug, which indeed was an edge-case bug, and I've > > >fixed it. Update to PSCBS v0.35.0: > > > > > >source(" >http://aroma-project.org/hbLite.R >"); > > >hbInstall("PSCBS"); > > > > > >I've verified that it works with your data. Below is the script I used. > > > > > >stopifnot(packageVersion("PSCBS") >= "0.35.0"); > > > > > >library("PSCBS"); > > > > > >library("R.devices"); > > >devOptions("png", width=1024); > > >setOption("devEval/args/force", FALSE); > > > > > >sampleName <- "ingrid"; > > > > > >## Load the data from a previously saved segmentation object. > > >pathname <- sprintf("%s.Rbin", sampleName); > > >fit <- loadObject(pathname); > > >data <- getLocusData(fit); > > >rm(fit); > > > > > > > > >## Locate centromeres and other large gaps in data > > >gaps <- findLargeGaps(data, minLength=1e6); > > >knownSegments <- gapsToSegments(gaps); > > >## => dim(knownSegments) == c(64, 4) > > > > > >## 'betaT' is already normalized using TumorBoost => tbn=FALSE. > > >fit <- segmentByNonPairedPSCBS(data, knownSegments=knownSegments, > > >avgDH="median", seed=0xBEEF, verbose=-10); > > >printf("Number of segments: %d\n", nbrOfSegments(fit)); > > >## Number of segments: 568 > > > > > >toPNG(sampleName, tags="tracks,avgDH=median", aspectRatio=0.6, { > > >plotTracks(fit); > > >}); > > > > > > > > >## Estimate global background level in [0,1] (due to normal > > >contamination and more) > > >kappa <- estimateKappa(fit, verbose=-10); > > >printf("Kappa: %g\n", kappa); > > >## Kappa: 0.374814 > > > > > > > > >## Call allelic balance > > >## (a) Estimate DH threshold for calling AB > > >deltaAB <- estimateDeltaAB(fit, kappa=kappa); # If skipped, will be > > >done internally > > >printf("Delta_AB: %g\n", deltaAB); > > >## Delta_AB: 0.146787 > > > > > >## (b) Call AB based on bootstrapped segment DH levels > > >fit <- callAB(fit, delta=deltaAB, verbose=-100); > > >printf("Number of AB segments: %d\n", sum(getSegments(fit)$abCall, >na.rm=TRUE)); > > >## Number of AB segments: 193 > > > > > >toPNG(sampleName, tags="tracks,avgDH=median,AB", aspectRatio=0.6, { > > >plotTracks(fit); > > >}); > > > > > > > > >## Call loss of heterozygosity (LOH) > > >## (a) Estimate C_1 threshold for calling LOH > > >deltaLOH <- estimateDeltaLOH(fit); # If skipped, will be done internally > > >printf("Delta_LOH: %g\n", deltaLOH); > > >## Delta_LOH: 0.566545 > > > > > >## (b) Call LOH based on bootstrapped segment C_1 levels > > >fit <- callLOH(fit, delta=deltaLOH, verbose=-10); > > >printf("Number of LOH segments: %d\n", sum(getSegments(fit)$lohCall, > > >na.rm=TRUE)); > > >## Number of LOH segments: 149 > > > > > >toPNG(sampleName, tags="tracks,avgDH=median,AB+LOH", aspectRatio=0.6, { > > >plotTracks(fit); > > >}); > > > > > > > > >## Call NTCN > > >## (a) Estimate the threshold for calling neutral TCN segments > > >## By shrinking 'scale', more segments will be non-NTCN. > > >deltaCN <- estimateDeltaCN(fit, scale=1.0, kappa=kappa); > > >printf("Delta_CN: %g\n", deltaCN); > > >## Delta_CN: 0.312593 > > > > > >## (b) Call NTCN based on bootstrapped segment TCN levels > > >fit <- callNTCN(fit, delta=deltaCN, verbose=-20); > > >printf("Number of NTCN segments: %d\n", sum(getSegments(fit)$ntcnCall, > > >na.rm=TRUE)); > > >## Number of NTCN segments: 303 > > > > > >toPNG(sampleName, tags="tracks,avgDH=median,AB+LOH+NTCN", >aspectRatio=0.6, { > > >plotTracks(fit); > > >}); > > > > > >/Henrik > > > > > > > > >On Mon, Apr 22, 2013 at 6:00 PM, Henrik Bengtsson > > > > > >You received this message because you are subscribed to the Google Group >aroma-affymetrix. >You can >post via email <mailto:aroma-affymetrix@googlegroups.com>. >To unsubscribe from this group, >send <mailto:aroma-affymetrix+unsubscr...@googlegroups.com> an empty >message. >For more options, >visit <http://groups.google.com/group/aroma-affymetrix/topics> this group. > >-- >-- >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/ <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.
