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
On Mon, Feb 18, 2013 at 4:38 AM, Oscar Rueda <oscar.ru...@cruk.cam.ac.uk> wrote:
> Thanks Henrik for all your help
>
>
> Yes, I'm using 0.32.5 and I called
>
>> fit <- segmentByPairedPSCBS(X, knownSegments=knownSegments, seed=4874,
>>verbose=-9, tbn=FALSE, avgDH="median")
>
> Cheers,
> Oscar
>
> On 14/02/2013 10:04, "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() (Was: Re: [aroma.affymetrix] Digest for
>>aroma-affymetrix@googlegroups.com - 2 Messages in 1 Topic)
>><#group_thread_0> [1 Update]
>>
>>
>> PSCBS and callNTCN() (Was: Re: [aroma.affymetrix] Digest for
>>aroma-affymetrix@googlegroups.com - 2 Messages in 1 Topic)
>><http://groups.google.com/group/aroma-affymetrix/t/5a9fab682eebac1e>
>>
>>Henrik Bengtsson <henrik.bengts...@aroma-project.org>
>> Feb 13 01:48PM -0800
>>
>>
>>
>>
>>Hi Oscar,
>>
>>
>>
>>
>>
>>I'm rather swamped right now, but just a few quick comments. Calling
>>
>>
>>so called copy-neutral segments via callNTCN() is a rather and
>>
>>
>>somewhat experimental method (I should state in the vignette). In
>>
>>
>>other words, it has not undergone equally much rigorous validation as
>>
>>
>>the other callers. Having said this, this caller is very of interest
>>
>>
>>to us (and also right now through several projects), so any feedback
>>
>>
>>is very useful. I will try to have a look at your data that you've
>>
>>
>>sent me offline, but I cannot promise when. From the error message,
>>
>>
>>it could simply be a bug rather than a problem with the
>>
>>
>>estimators/data.
>>
>>
>>
>>
>>
>>About the avgDH="median". Exactly what commands did you call? Also,
>>
>>
>>you're using PSCBS v0.32.5, correct?
>>
>>
>>
>>
>>
>>/Henrik
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>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.
>
>
--
--
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.
