[BioC] crlmm with Illumina arrays

[Lavinia Gordon]

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