Re: [aroma.affymetrix] PSCBS and callNTCN()

[Oscar Rueda]

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-af...@googlegroups.com"
<aroma-af...@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-af...@googlegroups.com - 2 Messages in 1 Topic)
><#group_thread_0> [1 Update]
>
>
> PSCBS and callNTCN() (Was: Re: [aroma.affymetrix] Digest for
>aroma-af...@googlegroups.com - 2 Messages in 1 Topic)
><http://groups.google.com/group/aroma-affymetrix/t/5a9fab682eebac1e>
>
>Henrik Bengtsson <henrik.b...@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-af...@googlegroups.com>.
>To unsubscribe from this group,
>send <mailto:aroma-affymetr...@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-af...@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-affymetr...@googlegroups.com.
>For more options, visit https://groups.google.com/groups/opt_out.
>
>
>

[Henrik Bengtsson]

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