[BioC] crlmm with Illumina arrays

[Benilton Carvalho]

Hi Lavinia,

how many samples are you using in this particular example?

thanks,

benilton

On Wed, Feb 17, 2010 at 11:22 PM, Lavinia Gordon
<lavinia.gordon at mcri.edu.au> wrote:
> To authors/users of crlmm (and probably VanillaICE)
>
> I am using crlmm with Illumina arrays (human610quadv1b)
> I have run into a few problems and would appreciate any tips on how to
> troubleshoot.
> The issues may be related to the quality of the arrays, however they have
> all run through the standard Illumina software without problems.
>
> I am following a (slightly) modified version of the script
> illumina_copynumber.Rnw,
> finishing with:
> crlmmWrapper(sampleSheet=samplesheet,
>         arrayNames=arrayNames,
>         arrayInfoColNames=list(barcode="SentrixBarcode_A",
> position="SentrixPosition_A"),
>         saveDate=TRUE,
>         cdfName=cdfName,
>         load.it=FALSE,
>         save.it=FALSE,
>         intensityFile=file.path(outdir, "normalizedIntensities.rda"),
>         crlmmFile=file.path(outdir, "snpsetObject.rda"),
>         rgFile=file.path(outdir, "rgFile.rda"))
>
> If I run:
> update(filename, adj.bias=TRUE)
> I get this error:
> Processing  OUTDIR/crlmmSetList_1.rda ...
> ----------------------------------------------------------------------------
> -        Estimating copy number for chromosome 1
> ----------------------------------------------------------------------------
> Error in computeCopynumber(object, CHR = CHR, ...) :
>  formal argument "CHR" matched by multiple actual arguments
> (not just for chromosome 1)
>
> So instead I have been using 'computeCopyNumber':
> crlmmSetList <- computeCopynumber(crlmmSetList, CHR, bias.adj=TRUE,
> SNRmin=5, cdfName, batch=scanbatch)
> save(crlmmSetList, file=file.path(outdir, paste("crlmmSetList_", CHR,
> ".rda", sep="")))
> This runs for all bar chromosome 23 and 24 (which is probably due to the
> fact that these .rda files are tiny, ~1.9kb)
> [1] "OUTDIR/crlmmSetList_23.rda"
> Fitting model for copy number estimation...
> Using 50 df for inverse chi squares.
> Estimating gender
> Error in kmeans(XMedian, c(min(XMedian[SNR > SNRmin]), max(XMedian[SNR >  :
>  initial centers are not distinct
>> traceback()
> 5: stop("initial centers are not distinct")
> 4: kmeans(XMedian, c(min(XMedian[SNR > SNRmin]), max(XMedian[SNR >
>       SNRmin])))
> 3: instantiateObjects(calls = calls, conf = conf, NP = NP, plate = plate,
>       envir = envir, chrom = chrom, A = A, B = B, gender = gender,
>       SNR = SNR, SNRmin = SNRmin, pkgname = cdfName)
> 2: .computeCopynumber(chrom = CHR, A = A(ABset), B = B(ABset), calls =
> calls(snpset),
>       conf = confs(snpset), NP = A(NPset), plate = batch, envir = envir,
>       SNR = ABset$SNR, bias.adj = FALSE, SNRmin = SNRmin, cdfName = cdfName,
>       ...)
> 1: computeCopynumber(crlmmSetList, CHR, bias.adj = TRUE, SNRmin = 5,
>       cdfName, batch = scanbatch)
>
> Then following a slightly modified version of the script copynumber.Rnw,
> beginning from the section:
> CHR <- 2
> if(!exists("crlmmSetList")) load(file.path(outdir, paste("crlmmSetList_",
> CHR, ".rda", sep="")))
> show(crlmmSetList)
> ....
> to
>  if(!exists("hmmPredictions")){
> + hmmPredictions <- viterbi(emission=emission.cn,
> +   initialStateProbs=log(initialP),
> +   tau=tau[, "transitionPr"],
> +   arm=tau[, "arm"],
> +   normalIndex=3,
> +   normal2altered=0.1,
> +   altered2normal=1,
> +   altered2altered=0.01)
> + }
>
> I get:
>
> Error: NA/NaN/Inf in foreign function call (arg 1)
>> summary(is.na(emission.cn))
>   Mode   FALSE    TRUE    NA's
> logical 8806640 3473576       0
>> head(tau)
>     chromosome position arm transitionPr
> [1,]          2     8856   0    0.9998882
> [2,]          2    14445   0    0.9998708
> [3,]          2    20906   0    0.9999908
> [4,]          2    21366   0    0.9999915
> [5,]          2    21791   0    0.9999756
> [6,]          2    23012   0    0.9999245
>> traceback()
> 2: .C("viterbi", tmp[[1]], tmp[[2]], tmp[[3]], tmp[[4]], tmp[[5]],
>       tmp[[6]], tmp[[7]], tmp[[8]], tmp[[9]], tmp[[10]], tmp[[11]],
>       tmp[[12]], tmp[[13]])
> 1: viterbi(emission = emission.cn, initialStateProbs = log(initialP),
>       tau = tau[, "transitionPr"], arm = tau[, "arm"], normalIndex = 3,
>       normal2altered = 0.1, altered2normal = 1, altered2altered = 0.01)
>> sessionInfo()
> R version 2.10.1 (2009-12-14)
> x86_64-unknown-linux-gnu
>
> 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=C              LC_MESSAGES=en_US.UTF-8
>  [7] LC_PAPER=en_US.UTF-8       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] human610quadv1bCrlmm_1.0.1 RColorBrewer_1.0-2
> [3] SNPchip_1.10.0             oligoClasses_1.8.0
> [5] VanillaICE_1.8.0           crlmm_1.4.1
> [7] Biobase_2.6.1
>
> loaded via a namespace (and not attached):
>  [1] affyio_1.14.0        annotate_1.24.1      AnnotationDbi_1.8.1
>  [4] Biostrings_2.14.10   DBI_0.2-5            ellipse_0.3-5
>  [7] genefilter_1.28.2    IRanges_1.4.9        mvtnorm_0.9-8
> [10] preprocessCore_1.8.0 RSQLite_0.8-1        splines_2.10.1
> [13] survival_2.35-7      tools_2.10.1         xtable_1.5-6
>>
>
> Any advice greatly appreciated,
> with regards
>
> Lavinia Gordon
> Bioinformatics officer
> Murdoch Childrens Research Institute
> The Royal Children s Hospital, Flemington Road, Parkville, Victoria 3052,
> Australia
>
> _______________________________________________
> Bioconductor mailing list
> Bioconductor at stat.math.ethz.ch
> https://stat.ethz.ch/mailman/listinfo/bioconductor
> Search the archives:
> http://news.gmane.org/gmane.science.biology.informatics.conductor
>