I used getSegmentReadCountsFromBAM() to get a GRanges object from a list of tumor-normal bam files Then I used referencecn.mops() which triggered the following (french) error :

NOTE: The default parameters are adjusted for "tumor-vs-normal"!

Normalizing...

Starting local modeling, please be patient...  

Reference sequence:  1

Starting segmentation algorithm...

Using "DNAcopy" for segmentation.

Erreur dans nperm * alpha : argument non numérique pour un opérateur binaire

> sessionInfo()

R version 3.1.2 (2014-10-31)

Platform: x86_64-redhat-linux-gnu (64-bit)

locale:

 [1] LC_CTYPE=fr_FR.UTF-8       LC_NUMERIC=C              

 [3] LC_TIME=fr_FR.UTF-8        LC_COLLATE=fr_FR.UTF-8    

 [5] LC_MONETARY=fr_FR.UTF-8    LC_MESSAGES=fr_FR.UTF-8   

 [7] LC_PAPER=fr_FR.UTF-8       LC_NAME=C                 

 [9] LC_ADDRESS=C               LC_TELEPHONE=C            

[11] LC_MEASUREMENT=fr_FR.UTF-8 LC_IDENTIFICATION=C       

attached base packages:

[1] stats4    parallel  stats     graphics  grDevices utils     datasets 

[8] methods   base     

other attached packages:

 [1] DNAcopy_1.40.0       ggplot2_2.0.0        mclust_5.1          

 [4] cit.supervised_1.0.0 qvalue_1.43.0        actuar_1.1-10       

 [7] limma_3.22.7         survival_2.38-3      cit.utils_1.0.0     

[10] MASS_7.3-45          rgl_0.95.1441        cn.mops_1.12.0      

[13] GenomicRanges_1.18.4 GenomeInfoDb_1.2.5   IRanges_2.0.1       

[16] S4Vectors_0.4.0      Biobase_2.26.0       BiocGenerics_0.12.1 

[19] BiocInstaller_1.16.5

loaded via a namespace (and not attached):

 [1] Biostrings_2.34.1 bitops_1.0-6      colorspace_1.2-6  grid_3.1.2       

 [5] gtable_0.1.2      munsell_0.4.2     plyr_1.8.3        Rcpp_0.12.3      

 [9] Rsamtools_1.18.3  scales_0.3.0      snow_0.4-1        splines_3.1.2    

[13] tools_3.1.2       XVector_0.6.0     zlibbioc_1.12.0  

Dear Remy,

This is probably due to NAs in the read count matrix. Could you please check the distribution of read counts in the read counts matrix. You must have programmed something like:

X <- getSegmentReadCountsFromBAM(myBamFiles, ...)

Please check:

colSums(as.matrix(values(X)))
summary(as.matrix(values(X)))

If you do not see any obvious problems, such as samples without any reads, could you please send me the read count matrix via email: klambauer@bioinf.jku.at

You can export the read count matrix as csv or Rds:

write.csv(as.data.frame(X), file="myReadCounts.csv")

or:

saveRDS(X,file="myReadCounts.Rds")

Regards,

Günter

Dear Günther,

I managed to solve my problem. So, I happen to load the "scales" package from the CRAN (I'm not even sure why, but I do by default) and this package happens to have a function called alpha. Now at line 462 of the referencecn.mops, the function tests wether alpha is set and if not sets it to 0.01 ( if (!exists("alpha")){alpha <- 0.01} ). Evidently, alpha already existed in my session and is a function.

So setting alpha to 0.01 before running the analysis does the trick in my case. 

Now concerning the analysis,

1. I was just trying to make it work but I indeed run it on all the chromosomes by adding all the chromosome names to the refSeqName parameter.

2. I'm not sure what you mean by classical but these are pairs of tumor/normal samples. 

3. These are actually pure tumor samples (Xenograft with a strong process) but I'd love to know how to take purity into account for other analysis.

4. I've got no idea what the ploidy of the tumor samples is, but I'd love to have an estimate from sequencing and has for point 3, I might get an estimate in other projects (or later on this one).

Did I miss something in the vignette or in the help of the function and I should know how to input all the information I've got in the cn.mops process?

Thanks for your input anyway.

Best,

Rémy

Thanks for your reply Günther.

I used getReadCountsFromBAM instead of getSegmentReadCountsFromBAM.

I checked the two matrices and all samples have more than 4M counts.

I removed every line in the matrices that had a zero value for all samples, about 4K lines in 50K total and I still got the same error.

I'll send you the files.

Remy

Dear Remy,

I could not reproduce your problems. referencecn.mops ran through cleanly giving reasonable results:

library(cn.mops)

load("normalbamDataRanges.RData")
load("tumorbamDataRanges.RData")

rr <- referencecn.mops(controls=normalbamDataRanges, cases=tumorbamDataRanges)
rr <- calcIntegerCopyNumbers(rr)

segplot(rr,sampleIdx=6)

write.csv(as.data.frame(cnvs(rr)),file="cnvs.csv")

Could you give me the code that gave the error?

A few more comments/questions:

1. Is this only chromosome 1 data? If "yes", I suggest you use chr1-22 together to get optimal results
2. Is this a classical tumor vs normal study?
3. Are these tumor cell lines or tumor tissue? Tumor tissue usually has impurity that needs to be taken into taken account with cn.mops parameters... Do you have purity estimates?
4. Non-diploidy of tumor also has to be taken into account with cn.mops parameters.

Let me know if you need further help with the analysis!

Regards,

Günter

Dear Remy,

Thanks for the update - this is indeed a very "unfortunate" problem with your code.

As for your analysis: I will contact you with details via Email!

Regards,

Günter