search

JEFworks-Lab / HoneyBADGER

HMM-integrated Bayesian approach for detecting CNV and LOH events from single-cell RNA-seq data
http://jef.works/HoneyBADGER/
GNU General Public License v3.0
96 stars 31 forks source link

Error in calcAlleleCnvBoundaries #26

Open JLTrincado opened 5 years ago

JLTrincado commented 5 years ago

Dear Jean,

first, congratulations for this method and the paper.

I'm testing it on a set of Smart-seq2 data for which we have also WGS information available. I extracted the counts for the reference and allele count and follow your tutorial for the downstream steps. I'm able to get the allele profile, but when I try to run calcAlleleCnvBoundaries I get the following error:

hb_04$calcAlleleCnvBoundaries(init=TRUE, verbose=TRUE) ## HMM ignore previously identified CNVs ... iterative HMM ... max vote:3 SNPS AFFECTED BY DELETION/LOH: chr11:36463422:36463422 chr11:68588015:68588015 chr11:107429423:107429423 chr11:108917891:108917891 chr11:114189367:114189367 chr11:134149641:134149641 chr12:909319:909319 chr12:10612961:10612961 chr12:27800374:27800374 chr12:54231632:54231632 chr12:54234519:54234519 chr12:56270842:56270842 chr12:96194399:96194399 chr12:109852372:109852372 chr12:120701764:120701764 chr12:131921676:131921676 chr13:21495681:21495681 chr13:23355937:23355937 chr13:27426307:27426307 chr13:32586535:32586535 chr13:39038397:39038397 chr13:52687852:52687852 chr14:74302884:74302884 chr14:92513948:92513948 chr14:95158667:95158667 chr14:103590772:103590772 chr15:50370390:50370390 chr15:63616479:63616479 chr15:64365021:64365021 chr15:64365340:64365340 chr16:11691753:11691753 chr16:24919799:24919799 chr16:71745491:71745491 chr17:5434694:5434694 chr17:7311357:7311357 chr17:28906269:28906269 chr17:32479946:32479946 chr17:43169904:43169904 chr17:50360705:50360705 chr17:63928898:63928898Assessing posterior probability of CNV in region ... with 40 snps ... within 38 genes ... converting to multi-dimensional arrays ... Error in n.sc.array[i, s, ] <- n.sc[snpst[s], ] : incorrect number of subscripts

I don't get why this error with the dimensions is happening. If I run the same but saving it on a separate object it seems to work, but the method runs differently:

potentialCnvs <- calcAlleleCnvBoundaries(allele.mats$r.maf,allele.mats$n.sc,allele.mats$l.maf,allele.mats$n.bulk, allele.mats$snps, geneFactor, verbose = TRUE) Running single round of HMM with tree cut height 3... max vote:3 [1] "Identified 20 potential CNVs"

Thanks a lot for your time. Best. Juan L.

JLTrincado commented 5 years ago

I attach my session info, in case it helps:

sessionInfo() R version 3.6.0 (2019-04-26) Platform: x86_64-pc-linux-gnu (64-bit) Running under: Ubuntu 18.04.2 LTS

Matrix products: default BLAS: /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.7.1 LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.7.1

locale: [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C LC_TIME=es_ES.UTF-8 LC_COLLATE=en_US.UTF-8 LC_MONETARY=es_ES.UTF-8
[6] LC_MESSAGES=en_US.UTF-8 LC_PAPER=es_ES.UTF-8 LC_NAME=es_ES.UTF-8 LC_ADDRESS=es_ES.UTF-8 LC_TELEPHONE=es_ES.UTF-8
[11] LC_MEASUREMENT=es_ES.UTF-8 LC_IDENTIFICATION=es_ES.UTF-8

attached base packages: [1] parallel stats4 stats graphics grDevices utils datasets methods base

other attached packages: [1] TxDb.Hsapiens.UCSC.hg19.knownGene_3.2.2 gridExtra_2.3 reshape2_1.4.3 TxDb.Hsapiens.UCSC.hg38.knownGene_3.4.6 [5] GenomicFeatures_1.36.0 AnnotationDbi_1.42.1 VariantAnnotation_1.30.0 Rsamtools_2.0.0
[9] Biostrings_2.52.0 XVector_0.23.2 SummarizedExperiment_1.14.0 DelayedArray_0.10.0
[13] BiocParallel_1.18.0 matrixStats_0.54.0 Biobase_2.40.0 GenomicRanges_1.35.2
[17] GenomeInfoDb_1.19.3 IRanges_2.17.6 S4Vectors_0.21.24 BiocGenerics_0.30.0
[21] HoneyBADGER_0.1 ggsignif_0.5.0 xlsx_0.6.1 ggplot2_3.1.1
[25] reshape_0.8.8 plyr_1.8.4

loaded via a namespace (and not attached): [1] httr_1.4.0 bit64_0.9-7 assertthat_0.2.1 BiocManager_1.30.4 HiddenMarkov_1.8-11 blob_1.1.1
[7] xlsxjars_0.6.1 BSgenome_1.52.0 GenomeInfoDbData_1.2.1 yaml_2.2.0 progress_1.2.0 pillar_1.3.1
[13] RSQLite_2.1.1 lattice_0.20-38 glue_1.3.1 digest_0.6.18 colorspace_1.4-1 htmltools_0.3.6
[19] Matrix_1.2-17 XML_3.98-1.19 pkgconfig_2.0.2 biomaRt_2.36.1 zlibbioc_1.29.0 purrr_0.3.2
[25] scales_1.0.0 rjags_4-8 tibble_2.1.1 withr_2.1.2 lazyeval_0.2.2 magrittr_1.5
[31] crayon_1.3.4 memoise_1.1.0 evaluate_0.13 tools_3.6.0 prettyunits_1.0.2 hms_0.4.2
[37] stringr_1.4.0 munsell_0.5.0 compiler_3.6.0 caTools_1.17.1.2 rlang_0.3.4 grid_3.6.0
[43] RCurl_1.95-4.12 rstudioapi_0.10 bitops_1.0-6 labeling_0.3 rmarkdown_1.12 gtable_0.3.0
[49] DBI_1.0.0 R6_2.4.0 GenomicAlignments_1.20.0 rtracklayer_1.44.0 knitr_1.22 dplyr_0.8.0.1
[55] bit_1.1-14 rJava_0.9-11 stringi_1.4.3 Rcpp_1.0.1 tidyselect_0.2.5 xfun_0.6
[61] coda_0.19-2

JLTrincado commented 5 years ago

Dear Jean,

sorry for insisting. Did you have time for taking a look on this?

Thanks again for all your work. Best, Juan L.

JEFworks commented 5 years ago

Hi Juan,

Thanks for your patience. This is on my task list and I will get to it as soon as possible.

My general sense that this error is related to a previously closed issue: https://github.com/JEFworks/HoneyBADGER/issues/13 There were a number of bugs like these that were fixed in the individual functions like calcAlleleCnvBoundaries() but were unfortunately not properly ported over to their object analogues ie. hb$calcAlleleCnvBoundaries