[R] FW: cgh package

[Tom Price]

Dear Celine,

Please excuse the delay in responding over the holiday period. Here
are my answer to your questions. If you are having further
difficulties or these answers don't seem right do let me know.

(1) Yes, I recommend splitting your data by chromosome. If you perform
a genomewide analysis, I recommend that you adjust the for testing
multiple chromosomes using the Bonferroni method. That is, for a
genomewide target p-value of alpha=0.05 you should employ a
chromosomewide p-value of alpha=0.05/14=0.0036.

(2) In your code, you should use

r <- sw.rob(x)

not

r <- sw.rob(intens[,8])

I should really update the documentation to make this clear.

The high scoring islands are contained in the results of sw(). The
blue line should start at sw$start[1] and have a length of
sw$length[1]. If this is not then there is probably some descrepency
between the output of sw() and what you are plotting. If you are still
having problems with this please email me with more details.

The sw$score is sorted with the highest scoring island first. The
start and lengths of the island given in sw$start and sw$length follow
the same order, as do the results of the permutation test (which are
p-values). As you point out the documentation does not mention this.

Tom



>
> Hi,
>
> I would like some extra information on the 'cgh' package in R. I noticed
> that there isn't much activity regarding this package on the R and BioC
> mailing list (I googled it).
>
>
>
> I started using this package and I have few questions:
>
> 1/ As I have a custom tiling like array @8um features resolution (affy), I
> have a lot of probes to work with. I'm assuming it is correct to split my
> data by chromosomes (I work on Plasmodium): 14 chromosomes, no sexual
> chromosome, so nothing really expected in terms of polysomy.
>
>
>
> 2/ For each chromosome, I'm computing the robustness, the threshold, the sw
> and the sw.perm.test before to plot with sign=-1 and a second plot sign=+1
>
> -          Do I need to change the threshold's sign when I plot?  (if the
> threshold is -1, I plot sign=-1, but if I do sign=+1 with the same threshold
> -1, how wrong is that?)
>
> -          How do you export the high scoring islands? If I inspect the sw
> object created by sw(), the highest score doesn't seem to correspond to the
> blue line generated by the sw.plot as "highest scoring island".
>
> -          In the sw object, are the values of score, start, length, sorted
> for the highest score? They seem to be… but I couldn't find any
> documentation about it
>
> -          Is the sw.perm also sorted for the highest score? Or is it for
> the highest p value?
>
>
>
> What I really want is to be able to write.table() for the islands location,
> and the best score, ideally, the ones corresponding to the island(s) plotted
> by the sw.plot().
>
>
>
> Here is the exemple code of what I used:
>
>
>
> x <- sw.threshold(intens[,8], function(x) median(x) + .2 * mad(x), -1)
>
> r <- sw.rob(intens[,8])
>
> sw <- sw(x,max.nIslands = NULL, trace = F)
>
> sw.perm <- sw.perm.test(x, max.nIslands = NULL, nIter= 1e4)
>
> sw.plot(intens[,8], location = coords[,2],threshold.func = function(x)
> median(x) + .2 * mad(x), sign = -1, expected = NULL, legend =
> TRUE,rob=r,main=paste("chrom 9 loss p=",sw.perm[1]))
>
> sessionInfo()
>
> R version 2.6.0 (2007-10-03)
>
> ia64-unknown-linux-gnu
>
>
>
> locale:
>
> C
>
>
>
> attached base packages:
>
> [1] stats     graphics  grDevices utils     datasets  methods   base
>
>
>
> other attached packages:
>
> [1] cgh_1.0-2
>
>
>
> loaded via a namespace (and not attached):
>
> [1] rcompgen_0.1-17
>
>
>
> Thank you very much for your help,
>
> Best wishes
>
> celine
>
>
>
> **************************************
>
> Celine Carret, PhD
>
> Pathogen Microarrays
>
> Wellcome Trust Sanger Institute
>
> Hinxton, Cambridge
>
> CB10 1SA, UK
>
> tel.+44-1223494940
>
> fax.+44-1223494919
>
> **************************************
>
>
>  -- The Wellcome Trust Sanger Institute is operated by Genome Research
> Limited, a charity registered in England with number 1021457 and a company
> registered in England with number 2742969, whose registered office is 215
> Euston Road, London, NW1 2BE.
>