Hi Ingrid.
[cleaning up by dropping everything but the
bootstrapTCNandDHByRegion() error; will reply to the other stuff
separately.]
On Sun, Apr 21, 2013 at 8:33 PM, Ingrid Lönnstedt
<ingrid.lonnst...@gmail.com> wrote:
>
> Hi Henrik,
> Thank you for your answers! You do a great job, indeed! I am now back to
> this. See a couple of remaining issues below.
> /Ingrid
>
>
>
> On Wednesday, 27 March 2013 18:16:07 UTC+11, Henrik Bengtsson wrote:
>>
>> On Tue, Mar 26, 2013 at 10:03 PM, Ingrid Lönnstedt
>> <ingrid.l...@gmail.com> wrote:
>> >
>> > Thank you Henrik!
>> > OK, I followed your instructions, and the same problem (the error in
>> > callAB()) remains, or rather, it now occurs already in sample one. Also,
[snip]
> So, my main problem that remains is insanity check failures in callAB(). I
> have followed my sample 1 into bootstrapTCNandDHByRegion.PairedPSCBS with
> debug() and concluded that I have 592 segments, and the first failure occurs
> with the 233rd segment. I have printed the CBS info for segments 220 to 233,
> in case that would help, as well as the current error message and session
> Info. If there is a failure with segment 233, is there a reason not to
> report the callAB outcomes of the other segments? What does the sanity
> check failure tell me, and could I do any fix to get round it or at least
> get something for the other segments?
>
> Thanks,
> Ingrid
>
>
> Browse[3]> segs[220:233,]
>
> chromosome tcnId dhId tcnStart tcnEnd tcnNbrOfLoci tcnMean
> 220 9 11 1 73187382 125406672 38303 2.1977136
> 221 9 12 1 125406672 141091395 8913 2.1506696
> 222 NA NA NA NA NA 0 NA
> 223 10 1 1 72760 17932882 15166 2.0003566
> 224 10 2 1 17932882 17956182 2 4.5480683
> 225 10 3 1 17956182 18179544 18 1.9775509
> 226 10 4 1 18179544 18202780 3 0.6127961
> 227 10 5 1 18202780 24376470 4403 2.0286301
> 228 10 6 1 24376470 24378482 11 0.6265435
> 229 10 7 1 24378482 39076222 10124 2.0268394
> 230 10 8 1 39076223 42433539 0 NaN
> 231 10 9 1 42433540 86411230 28982 2.0375534
> 232 10 10 1 86411230 89029244 1697 2.0558821
> 233 10 11 1 89029244 89115851 21 2.5505440
> tcnNbrOfSNPs tcnNbrOfHets dhNbrOfLoci dhStart dhEnd dhMean
> 220 18928 5290 5290 73187382 125406672 0.07756713
> 221 4002 1101 1101 125406672 141091395 0.07943612
> 222 NA NA 0 NA NA NA
> 223 8456 2087 2087 72760 17932882 0.07487577
> 224 0 0 0 17932882 17956182 NA
> 225 4 1 1 17956182 18179544 0.04524338
> 226 0 0 0 18179544 18202780 NA
> 227 2212 543 543 18202780 24376470 0.07859021
> 228 0 0 0 24376470 24378482 NA
> 229 4971 1086 1086 24378482 39076222 0.07712847
> 230 0 0 0 39076223 42433539 NA
> 231 13983 3576 3576 42433540 86411230 0.06990635
> 232 858 192 192 86411230 89029244 0.06439847
> 233 4 1 1 89029244 89115851 0.36474615
> c1Mean c2Mean
> 220 1.0136216 1.184092
> 221 0.9899144 1.160755
> 222 NA NA
> 223 0.9252892 1.075067
> 224 NA NA
> 225 0.9440399 1.033511
> 226 NA NA
> 227 0.9345998 1.094030
> 228 NA NA
> 229 0.9352562 1.091583
> 230 NA NaN
> 231 0.9475577 1.089996
> 232 0.9617432 1.094139
> 233 0.8101214 1.740423
>
> deltaCN <- estimateDeltaCN(fit, scale=1) #0.2601609 for my sample
> fit = callAB(fit, delta = deltaCN, verbose = -10)
> List of 7
> $ nbrOfTCNs : int 21
> $ tcnNbrOfLoci : int 21
> $ mu : num 2.58
> $ tcnMean : num 2.55
> $ dMu : num 0.0248
> $ abs(dMu) : num 0.0248
> $ range(x[units]): num [1:2] 89029244 89114401
>
> Error in bootstrapTCNandDHByRegion.PairedPSCBS(fit, statsFcn = statsFcn,
> :
> INTERNAL ERROR: Incorrect TCN mean!
This is an internal error that the *implementation* of the PSCBS
algorithm is correct and handles all "edge cases"; apparently there is
still another case that was not considered. These "edge cases" are
for instance when two SNPs share the exact same location but the
internal CBS segmentation method wants to put a change point in
between them because that's what the data happens to say. In other
words, it's nothing you're doing wrong and probably nothing you should
try to work around either. (I'll think about this a bit more since it
happens once in a while and it would be nice to be able to keep going.
On the other hand, if it would be just a warning you would have never
told me and the error would never be fixed.)
First, update to the latest PSCBS v0.34.8 (on CRAN since yday) and
verify that the error remains. If it does, could you save the
segmentation result;
saveObject(fit, "problem.Rdata")
and make 'problem.Rdata' available to me for download somewhere? Then
I can troubleshoot.
/Henrik
>
>
> > sessionInfo()
> R version 2.15.2 Patched (2013-02-09 r61879)
> 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] stats graphics grDevices utils datasets methods base
>
> other attached packages:
> [1] DNAcopy_1.32.0 PSCBS_0.34.1 aroma.affymetrix_2.8.0
> [4] affxparser_1.30.2 aroma.apd_0.2.3 R.huge_0.4.1
> [7] aroma.light_1.28.0 aroma.core_2.8.0 matrixStats_0.7.0
> [10] R.rsp_0.8.2 R.devices_2.1.3 R.cache_0.6.5
> [13] R.filesets_2.0.0 R.utils_1.23.2 R.oo_1.13.0
> [16] R.methodsS3_1.4.2
>
> loaded via a namespace (and not attached):
> [1] calmate_0.10.0 digest_0.6.3 MASS_7.3-23
>
>
>>
>> > All the best,
>> > Ingrid
>> >
>> >
>> > On Saturday, 23 March 2013 05:03:50 UTC+11, Henrik Bengtsson wrote:
>> >>
>> >> Hi Ingrid.
>> >>
>> >> On Wed, Mar 20, 2013 at 4:45 PM, Ingrid Lönnstedt
>> >> <ingrid.l...@gmail.com> wrote:
>> >> >
>> >> > Hi!
>> >> >
>> >> > I have segmented 56 SNP 6 tumour samples with
>> >> > segmentByNonPairedPSCBS(),
>> >> > and like to call allelic balance segments with callAB() as required,
>> >> > I
>> >> > understand, before I can call neutral TCNs with callNTCN().
>> >>
[snip]
>> >>
>> >> source("http://aroma-project.org/hbLite.R";);
>> >> hbInstall("PSCBS");
>> >>
>> >>
>> >>
>> >> > For the first
>> >> > tumour sample, both callAB() and callNTCN() seem to work fine, but
>> >> > then on
>> >> > the second sample I get errors with callAB(), see below.
>> >>
>> >> See comment after the error.
>> >>
>> >> > The plotTracks()
>> >> > image looks pretty much the same between sample one and two, and I
>> >> > have no
>> >> > reason to believe that sample two has a major problem in the
>> >> > segmentation. 4
>> >> > of the 56 samples get very noisy segmentations, with lots of short
>> >> > regions
>> >> > scattered (has anybody seen that before?)
[snip]
>> >> > betaT = dataN[,'fracB',samp]
>> >> > ugp = getAromaUgpFile(dsNList$total)
>> >> > chromosome = ugp[,1,drop = T]
>> >> > x = ugp[, 2, drop = T]
>> >> > df = data.frame(chromosome = chromosome, x=x, CT=CT, betaT =
>> >> > betaT)
>> >> > df = dropSegmentationOutliers(df)
>> >> > #This gives warnings:
>> >> > #24: In DNAcopy::CNA(genomdat = yKK, chrom = chromosomeKK,
>> >> > ... :
>> >> > #array has repeated maploc positions
>> >> > gaps = findLargeGaps(df, minLength = 1e+06)
>> >> > knownSegments = gapsToSegments(gaps)
>> >> >
>> >> > fit = segmentByNonPairedPSCBS(df, knownSegments = knownSegments,
>> >> > avgDH="median",
>> >> > tauA = .33,
>> >> > seed =
>> >> > 48879, verbose = -10)
>> >> > #This gives warnings that there are multipe maploc positions
>> >> > fit = callAB(fit, verbose = -10)
>> >> >
>> >> > Calling segments of allelic balance from one-sided DH bootstrap
>> >> > confidence
>> >> > intervals...
>> >> > delta (offset adjusting for bias in DH): 0.211106073856354
>> >> > alpha (CI quantile; significance level): 0.05
>> >> > List of 7
>> >> > $ nbrOfDHs : int 70
>> >> > $ dhNbrOfLoci : int 70
>> >> > $ mu : num 0.183
>> >> > $ dhMean : num 0.203
>> >> > $ dMu : num -0.02
>> >> > $ abs(dMu) : num 0.02
>> >> > $ range(x[units]): num [1:2] 990417 2792801
>> >> > Error in bootstrapTCNandDHByRegion.PairedPSCBS(fit, statsFcn =
>> >> > statsFcn,
>> >> > :
>> >> > INTERNAL ERROR: Incorrect DH mean!
>> >> > Calling segments of allelic balance from one-sided DH bootstrap
>> >> > confidence intervals...done
>> >>
>> >> This is a internal sanity check throwing an error. As a start, I
>> >> recommend to retry with (i) a proper CT calculate, and (ii) after
>> >> updated PSCBS.
>> >>
>> >> Hope this helps
>> >>
>> >> Henrik
>> >>
>> >> >
>> >> >
>> >> > >
>> >> > > head(getSegments(fit))
>> >> > chromosome tcnId dhId tcnStart tcnEnd tcnNbrOfLoci tcnMean
>> >> > tcnNbrOfSNPs tcnNbrOfHets dhNbrOfLoci dhStart dhEnd
>> >> > 1 1 1 1 61736.0 811403.5 73 1.4735545
>> >> > 15 2 2 61736.0 811403.5
>> >> > 2 1 2 1 811403.5 2821394.0 527 0.9201369
>> >> > 183 70 70 811403.5 2821394.0
>> >> > 3 1 3 1 2821394.0 2827143.5 2 5.0273015
>> >> > 2 0 0 2821394.0 2827143.5
>> >> > 4 1 4 1 2827143.5 11694138.0 5531 0.9780667
>> >> > 2940 983 983 2827143.5 11694138.0
>> >> > 5 1 5 1 11694138.0 11709225.0 5 3.3949294
>> >> > 5 3 3 11694138.0 11709225.0
>> >> > 6 1 6 1 11709225.0 16923050.0 2918 1.0072282
>> >> > 1485 550 550 11709225.0 16923050.0
>> >> > dhMean c1Mean c2Mean
>> >> > 1 0.1695309 0.6118708 0.8616838
>> >> > 2 0.2033105 0.3665317 0.5536052
>> >> > 3 NA NA NA
>> >> > 4 0.2033807 0.3895734 0.5884933
>> >> > 5 0.1465555 1.4486918 1.9462376
>> >> > 6 0.2036166 0.4010699 0.6061583
>> >> >
>> >> >
>> >> > > traceback()
>> >> > 7: stop("INTERNAL ERROR: Incorrect DH mean!")
>> >> > 6: bootstrapTCNandDHByRegion.PairedPSCBS(fit, statsFcn = statsFcn,
>> >> > ..., verbose = less(verbose, 50))
>> >> > 5: bootstrapTCNandDHByRegion(fit, statsFcn = statsFcn, ..., verbose
>> >> > =
>> >> > less(verbose,
>> >> > 50))
>> >> > 4: callAllelicBalanceByDH.PairedPSCBS(fit, ...)
>> >> > 3: callAllelicBalanceByDH(fit, ...)
>> >> > 2: callAB.PairedPSCBS(fit, verbose = -10)
>> >> > 1: callAB(fit, verbose = -10)
>> >> >
>> >> >
>> >> >
>> >> >
>> >> > > sessionInfo()
>> >> >
>> >> > R version 2.15.2 Patched (2013-02-09 r61879)
>> >> > Platform: x86_64-unknown-linux-gnu (64-bit)
>> >> >
>> >> > locale:
>> >> > [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
>> >> > 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
>> >> > LC_PAPER=C
>> >> > LC_NAME=C
>> >> > [9] LC_ADDRESS=C LC_TELEPHONE=C
>> >> > LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
>> >> >
>> >> > attached base packages:
>> >> > [1] splines stats graphics grDevices utils datasets
>> >> > methods
>> >> > base
>> >> >
>> >> > other attached packages:
>> >> > [1] Hmisc_3.10-1 survival_2.37-4 DNAcopy_1.32.0
>> >> > PSCBS_0.30.0 aroma.affymetrix_2.8.0
>> >> > [6] affxparser_1.30.2 aroma.apd_0.2.3 R.huge_0.4.1
>> >> > aroma.light_1.28.0 aroma.core_2.8.0
>> >> > [11] matrixStats_0.6.2 R.rsp_0.8.2 R.devices_2.1.3
>> >> > R.cache_0.6.5 R.filesets_2.0.0
>> >> > [16] crlmm_1.16.9 preprocessCore_1.20.0
>> >> > oligoClasses_1.20.0
>> >> > BiocGenerics_0.4.0 R.utils_1.19.5
>> >> > [21] R.oo_1.11.7 R.methodsS3_1.4.2
>> >> >
>> >> > loaded via a namespace (and not attached):
>> >> > [1] affyio_1.26.0 annotate_1.36.0 AnnotationDbi_1.20.5
>> >> > Biobase_2.18.0 BiocInstaller_1.8.3 Biostrings_2.26.3
>> >> > [7] bit_1.1-9 cluster_1.14.3 codetools_0.2-8
>> >> > DBI_0.2-5 digest_0.6.3 ellipse_0.3-7
>> >> > [13] ff_2.2-10 foreach_1.4.0 genefilter_1.40.0
>> >> > GenomicRanges_1.10.7 grid_2.15.2 IRanges_1.16.6
>> >> > [19] iterators_1.0.6 lattice_0.20-13 Matrix_1.0-11
>> >> > mvtnorm_0.9-9994 parallel_2.15.2 RcppEigen_0.3.1.2.1
>> >> > [25] RSQLite_0.11.2 stats4_2.15.2 tools_2.15.2
>> >> > XML_3.95-0.1 xtable_1.7-1 zlibbioc_1.4.0
>> >> >
>> >> >
>> >> > >
>> >> >
>> >> >
>> >> >
>> >> > --
>> >> > --
>> >> > 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/
>> >> >
>> >> > ---
>> >> > 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.
>> >> >
>> >> >
>> >
>> > --
>> > --
>> > 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/
>> >
>> > ---
>> > 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.
>> >
>> >
--
--
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.
