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
On Sat, Feb 9, 2013 at 9:47 AM, Oscar Rueda <oscar.ru...@cruk.cam.ac.uk> wrote:
> 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.ru...@cruk.cam.ac.uk
>
>
>>Error in callAB (PSCBS)
>><http://groups.google.com/group/aroma-affymetrix/t/d066a381760d0543>
>>
>>Henrik Bengtsson <henrik.bengts...@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.bengts...@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.bengts...@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-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.
