Error in callAB (PSCBS)

[Oscar Rueda]

Dear all,

I'm trying to run PSCBS on 126 matched pairs, but only 60 of them have run
without errors. The last error I got was when running callAB:

> deltaAB <- estimateDeltaAB(fit, scale=1)
> fit <- callAB(fit, delta=deltaAB, verbose=-10)
Calling segments of allelic balance from one-sided DH bootstrap confidence
intervals...
delta (offset adjusting for bias in DH): 0.181206399202347
alpha (CI quantile; significance level): 0.05
Error: all(segMean + tol >= range[, 1], na.rm = TRUE) is not TRUE
Calling segments of allelic balance from one-sided DH bootstrap confidence
intervals...done

> traceback()
17: base::stop(error)
16: throw.error(ex)
15: throw(ex)
14: value[[3L]](cond)
13: tryCatchOne(expr, names, parentenv, handlers[[1L]])
12: tryCatchList(expr, classes, parentenv, handlers)
11: tryCatch({
fields <- dimnames(S)[[3]]
for (kk in seq(along = fields)) {
field <- fields[kk]
verbose && enter(verbose, sprintf("Field #%d ('%s') of %d",
kk, field, length(fields)))
Skk <- S[, , kk, drop = FALSE]
dim(Skk) <- dim(Skk)[-3]
stopifnot(is.matrix(Skk))
range <- Skk[, c(1, ncol(Skk)), drop = FALSE]
key <- sprintf("%sMean", field)
segMean <- segs[[key]]
verbose && printf(verbose, "mean: %g [%g,%g]\n", segMean,
range[1], range[2])
verbose && printf(verbose, "mean: %g, range: [%g,%g]\n",
segMean, range[1], range[2])
cfield <- sprintf("%sNbrOfLoci", ifelse(field == "tcn",
"tcn", "dh"))
counts <- segs[, cfield, drop = TRUE]
keep <- (counts > 1)
range <- range[keep, , drop = FALSE]
segMean <- segMean[keep]
stopifnot(all(range[, 2] + tol >= range[, 1], na.rm = TRUE))
stopifnot(all(segMean + tol >= range[, 1], na.rm = TRUE))
stopifnot(all(segMean - tol <= range[, 2], na.rm = TRUE))
verbose && exit(verbose)
}
}, error = function(ex) {
verbose && cat(verbose, "Tolerance (option
'PSCBS/sanityChecks/tolerance'): ",
tol)
verbose && print(verbose, segs)
throw(ex)
})
10: bootstrapTCNandDHByRegion.PairedPSCBS(fit, statsFcn = statsFcn,
..., verbose = less(verbose, 50))
9: bootstrapTCNandDHByRegion(fit, statsFcn = statsFcn, ..., verbose =
less(verbose,
50))
8: callAllelicBalanceByDH.PairedPSCBS(fit, ...)
7: callAllelicBalanceByDH(fit, ...)
6: callAB.PairedPSCBS(fit, delta = deltaAB)
5: callAB(fit, delta = deltaAB) at PSCSB_Manually.R#103
4: eval(expr, envir, enclos)
3: eval(ei, envir)
2: withVisible(eval(ei, envir))
1: source("PSCSB_Manually.R")

> sessionInfo()
R version 2.15.1 (2012-06-22)
Platform: x86_64-unknown-linux-gnu (64-bit)

locale:
[1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
[3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
[5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
[7] LC_PAPER=C LC_NAME=C
[9] LC_ADDRESS=C LC_TELEPHONE=C
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C

attached base packages:
[1] splines stats graphics grDevices utils datasets methods
[8] base

other attached packages:
[1] Hmisc_3.10-1 survival_2.36-14 DNAcopy_1.32.0
[4] PSCBS_0.30.0 aroma.affymetrix_2.6.0 affxparser_1.30.0
[7] aroma.apd_0.2.3 R.huge_0.4.1 aroma.light_1.28.0
[10] aroma.core_2.8.0 matrixStats_0.6.2 R.rsp_0.8.2
[13] R.devices_2.1.3 R.cache_0.6.5 R.filesets_1.9.0
[16] R.utils_1.19.3 R.oo_1.11.7 R.methodsS3_1.4.2

loaded via a namespace (and not attached):
[1] cluster_1.14.3 digest_0.6.0 grid_2.15.1 lattice_0.20-10
[5] tools_2.15.1




Has anyone encountered a similar problem?


Thanks,
Oscar






Oscar M. Rueda, PhD.
Postdoctoral Research Fellow, Caldas Lab, Breast Cancer Functional
Genomics.
Cancer Research UK Cambridge Research Institute.
Li Ka Shing Centre, Robinson Way.
Cambridge CB2 0RE
England

Please note that my email address will be changing to:
Oscar...@cruk.cam.ac.uk

[Henrik Bengtsson]

Hi,

sorry for the delay - I've been "off the grid" since two weeks.

On Thu, Jan 24, 2013 at 4:09 AM, Oscar Rueda <Oscar...@cruk.cam.ac.uk> wrote:
> Dear all,
>
> I'm trying to run PSCBS on 126 matched pairs, but only 60 of them have run
> without errors. The last error I got was when running callAB:
>
>> deltaAB <- estimateDeltaAB(fit, scale=1)
>> fit <- callAB(fit, delta=deltaAB, verbose=-10)
> Calling segments of allelic balance from one-sided DH bootstrap confidence
> intervals...
> delta (offset adjusting for bias in DH): 0.181206399202347
> alpha (CI quantile; significance level): 0.05
> Error: all(segMean + tol >= range[, 1], na.rm = TRUE) is not TRUE
> Calling segments of allelic balance from one-sided DH bootstrap confidence
> intervals...done

That is an internal "sanity" check part of a bootstrap estimation that
fails. The code has a lot of such assertion tests to validate that
intermediate and final estimates are sound and valid. It may be that
this one is a bit too conservative.

To troubleshoot, with a sample that fails, first confirm that you get
the same error by calling the internal bootstrap estimator explicitly:

fitB <- bootstrapTCNandDHByRegion(fit, verbose=-10);

It should, but if not, try another sample (aka 'fit'). Then, look for
the verbose output that looks something like:

Field #2 ('dh') of 4...
mean: 0.6645 [0.648152,0.112718]
mean: 0.1282 [0.648152,0.112718]
mean: 0.3038 [0.648152,0.112718]
mean: 0.6645, range: [0.648152,0.112718]
mean: 0.1282, range: [0.648152,0.112718]
mean: 0.3038, range: [0.648152,0.112718]
...
Tolerance (option 'PSCBS/sanityChecks/tolerance'): 5e-04

and some segment table followed by the error message. What do you get?

Also, when you find a sample that gives the error, save it

saveObject(fit, "faultySample.Rbin")

and make the file available to me for download (<20Mb can be emailed)
and I can have a look at it.

/Henrik

> --
> 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/

[Henrik Bengtsson]

Actually, before doing this, update to PSCBS 0.32.4;

source("http://aroma-project.org/hbLite.R");
hbInstall("PSCBS");

This should give a verbose output that instead looks like:

Field #2 ('dh') of 4...

mean=0.6645, range=[0.64765,0.68129], n=1047
mean=0.1282, range=[0.114076,0.143011], n=388
mean=0.3038, range=[0.290933,0.317302], n=664

...
Tolerance (option 'PSCBS/sanityChecks/tolerance'): 5e-04

which is more informative for the troubleshooting. Now what do you get?

/Henrik

[Oscar Rueda]

Hi Henrik,

Thanks for your reply.
I get

Statistical sanity checks (iff B >= 100)...
Field #1 ('tcn') of 4...
mean=1.8047, range=[1.79672,1.81258], n=15156
mean=1.7103, range=[1.68683,1.73451], n=1449
mean=1.8148, range=[1.80518,1.8243], n=10614
mean=1.9063, range=[1.88086,1.93251], n=1757
mean=1.8063, range=[1.78646,1.82722], n=2609
mean=1.9068, range=[1.89572,1.91696], n=10634
mean=1.968, range=[1.95805,1.97863], n=11580
mean=1.9316, range=[1.92209,1.94038], n=14313
mean=1.5654, range=[1.44976,1.69253], n=24
mean=2.13815, range=[2.10301,2.16987], n=1111
mean=2.2063, range=[2.15613,2.24568], n=769
mean=2.17853, range=[2.16259,2.19569], n=5849
mean=1.3989, range=[0.713311,2.00918], n=5
mean=NA, range=[NA,NA], n=0
mean=2.4012, range=[2.33236,2.46487], n=352
mean=2.9706, range=[2.91773,3.03072], n=849
mean=2.4297, range=[2.36293,2.51297], n=345
mean=2.7913, range=[2.73951,2.83993], n=932
mean=2.9512, range=[2.90223,3.0001], n=1206
mean=2.7529, range=[2.71359,2.79243], n=1568
mean=5.0963, range=[4.38476,5.74181], n=6
mean=2.8039, range=[2.74449,2.86477], n=649
mean=2.9969, range=[2.98574,3.0094], n=17880
mean=3.546, range=[3.49118,3.59413], n=1278
mean=3.8656, range=[3.8268,3.90099], n=4193
mean=3.8656, range=[3.8268,3.90099], n=4193
mean=2.7748, range=[2.59764,2.95961], n=63
mean=3.1181, range=[3.08514,3.15304], n=2665
mean=3.3677, range=[3.30253,3.43419], n=898
mean=3.198, range=[3.12122,3.2734], n=559
mean=2.7179, range=[2.59645,2.84069], n=189
mean=3.1643, range=[3.11939,3.20684], n=1812
mean=2.9993, range=[2.93235,3.06373], n=644
mean=2.808, range=[2.78317,2.83229], n=3725
mean=2.9411, range=[2.90161,2.98048], n=1505
mean=3.1943, range=[3.12907,3.26271], n=649
mean=2.9318, range=[2.90794,2.95572], n=4118
mean=3.188, range=[3.15639,3.21879], n=3282
mean=3.0313, range=[2.99859,3.06585], n=2459
mean=2.8684, range=[2.84352,2.89278], n=4497
mean=3.0151, range=[2.9895,3.04142], n=3718
mean=2.7903, range=[2.74472,2.83803], n=1027
mean=3.1104, range=[2.9115,3.2219], n=94
mean=5.2423, range=[4.4511,6.34827], n=5
mean=2.9283, range=[2.86237,2.99806], n=511
mean=3.151, range=[3.11551,3.18474], n=2364
mean=2.972, range=[2.95164,2.99418], n=6345
mean=NA, range=[NA,NA], n=0
mean=2.1594, range=[2.13923,2.18116], n=3169
mean=2.2456, range=[2.22776,2.26367], n=4018
mean=2.1282, range=[2.10191,2.1552], n=1900
mean=3.9688, range=[3.61876,4.3539], n=5
mean=2.1258, range=[2.06446,2.19542], n=255
mean=2.2482, range=[2.23571,2.26092], n=8513
mean=2.1435, range=[2.12764,2.15836], n=5278
mean=2.2946, range=[2.27967,2.30901], n=7239




and many more that contains some NAs

I have sent you another email with a link to download the R object with
the fit.

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...@cruk.cam.ac.uk

>Error in callAB (PSCBS)
><http://groups.google.com/group/aroma-affymetrix/t/d066a381760d0543>
>
>Henrik Bengtsson <henrik.b...@aroma-project.org>
> Feb 07 05:33PM -0800

>
>
>
>
>Hi,
>
>
>
>
>
>sorry for the delay - I've been "off the grid" since two weeks.
>
>
>
>
>

>> Error: all(segMean + tol >= range[, 1], na.rm = TRUE) is not TRUE
>
>
>> Calling segments of allelic balance from one-sided DH bootstrap
>>confidence
>
>
>> intervals...done
>
>
>
>
>
>That is an internal "sanity" check part of a bootstrap estimation that
>
>
>fails. The code has a lot of such assertion tests to validate that
>
>
>intermediate and final estimates are sound and valid. It may be that
>
>
>this one is a bit too conservative.
>
>
>
>
>
>To troubleshoot, with a sample that fails, first confirm that you get
>
>
>the same error by calling the internal bootstrap estimator explicitly:
>
>
>
>
>
>fitB <- bootstrapTCNandDHByRegion(fit, verbose=-10);
>
>
>
>
>
>It should, but if not, try another sample (aka 'fit'). Then, look for
>
>
>the verbose output that looks something like:
>
>
>
>
>
>Field #2 ('dh') of 4...
>
>
>mean: 0.6645 [0.648152,0.112718]
>
>
>mean: 0.1282 [0.648152,0.112718]
>
>
>mean: 0.3038 [0.648152,0.112718]
>
>
>mean: 0.6645, range: [0.648152,0.112718]
>
>
>mean: 0.1282, range: [0.648152,0.112718]
>
>
>mean: 0.3038, range: [0.648152,0.112718]
>
>
>...
>
>
>Tolerance (option 'PSCBS/sanityChecks/tolerance'): 5e-04
>
>
>
>
>
>and some segment table followed by the error message. What do you get?
>
>
>
>
>

>Also, when you find a sample that gives the error, save it
>
>
>
>
>
>saveObject(fit, "faultySample.Rbin")
>
>
>
>
>
>and make the file available to me for download (<20Mb can be emailed)
>
>
>and I can have a look at it.
>
>
>
>
>
>/Henrik
>
>
>
>
>
>
>
>
>
>
>
>

>Henrik Bengtsson <henrik.b...@aroma-project.org>
> Feb 07 06:03PM -0800

>
>
>
>
>On Thu, Feb 7, 2013 at 5:33 PM, Henrik Bengtsson
>
>

>> ...
>
>
>> Tolerance (option 'PSCBS/sanityChecks/tolerance'): 5e-04
>
>
>
>
>
>> and some segment table followed by the error message. What do you get?
>
>
>
>
>
>Actually, before doing this, update to PSCBS 0.32.4;
>
>
>
>
>
>source("
>http://aroma-project.org/hbLite.R
>");
>
>
>hbInstall("PSCBS");
>
>
>
>
>
>This should give a verbose output that instead looks like:
>
>
>
>
>
>Field #2 ('dh') of 4...
>
>
>mean=0.6645, range=[0.64765,0.68129], n=1047
>
>
>mean=0.1282, range=[0.114076,0.143011], n=388
>
>
>mean=0.3038, range=[0.290933,0.317302], n=664
>
>
>...
>
>
>Tolerance (option 'PSCBS/sanityChecks/tolerance'): 5e-04
>
>
>
>
>
>which is more informative for the troubleshooting. Now what do you get?
>
>
>
>
>
>/Henrik
>
>
>
>
>
>
>
>

> AromaUnitTotalCnBinarySet$byName- Error in cdf[[1]] : subscript out of
>bounds
><http://groups.google.com/group/aroma-affymetrix/t/413453191153b6aa>
>
>Henrik Bengtsson <henrik.b...@aroma-project.org>
> Feb 07 05:15PM -0800
>
>
>
>
>>> ds_CHI <- AromaUnitTotalCnBinarySet$byName(csR_CHI, tags=tags,
>
>
>>> chipType=chipType)
>
>
>> Error in cdf[[1]] : subscript out of bounds
>
>
>
>
>
>Did you really mean to use object 'csR_CHI' in that last call, or did you
>mean:
>
>
>
>
>
>ds_CHI <- AromaUnitTotalCnBinarySet$byName("47Chinese_CEL", tags=tags,
>
>
>chipType=chipType)
>
>
>
>
>
>?
>
>
>
>
>
>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]

Hi Oscar,

so you spotted a bug in the bootstrapping of the TCN mean-level
estimates. It only affected segments that had been called to be in
run-of-homozygosity (ROH) resulting in sampling from not all available
TCN loci. The effect was incorrect TCN mean quantile estimates for
those segments. For reasonable large segments this didn't matter
much, but for tiny segments (say < 5 loci) it made a small difference.
In your case the internal sanity checks caught this error in an ROH
segment of two loci.

The bug has been fixed in PSCBS v0.32.5. You can grab it as usual:

source("http://aroma-project.org/hbLite.R");
hbInstall("PSCBS");

I've verified that it now works with the sample data you've sent me.
Let me know if it solves it for all 126-60=66 samples where you had
problems.


BTW, It is experimental/still early days, but you may want to move to
calculate DH using the median estimator instead of the default sample
mean. For more details, see 'Less biased Decrease of Heterozygosity
(DH) estimates' in Section 'Experimental' of the 'Paired PSCBS'
vignette of PSCBS v0.32.5. Also, if you have saved 'fit' objects, you
don't have to rerun the segmentation (=finding the change points) to
achieve this; you can do fit <- updateMeans(fit, avgDH="median") and
then redo the calling.


Thanks for reporting.

Henrik