chr1swallace
commented
2 years ago thank you for the bug report. I can't reproduce this with my example data though. I did just fix a different bug to do with passing the sample size to susie when the susie doesn't converge on the first run. Could I ask you to please check with the latest coloc 5.2.1 and see if this somehow fixes your issue too?
Thanks.
https://chr1swallace.github.io
From: Alejandro Martinez @.> Sent: 21 November 2022 12:56 To: chr1swallace/coloc @.> Cc: Subscribed @.***> Subject: [chr1swallace/coloc] N argument missing when calling susieR::susie_rss() (Issue #107)
Hi,
I just looked at the wrapper to call susie_rss().
Even though you add n as an argument to susie_args list, , and this should be properly passed to the do.call function wrapping susie_rss, when the function runs the n argument is missing.
-
The data
str(coloc_data$gwas) List of 8 $ beta : Named num [1:4827] -0.0605 0.2421 -0.0605 0.0843 0.3397 ... ..- attr(, "names")= chr [1:4827] "rs12046188" "rs115541804" "rs12044651" "rs11162246" ... $ varbeta : Named num [1:4827] 0.00615 0.04435 0.00615 0.00792 0.11979 ... ..- attr(, "names")= chr [1:4827] "rs12046188" "rs115541804" "rs12044651" "rs11162246" ... $ N : num 2664 $ type : chr "quant" $ snp : chr [1:4827] "rs12046188" "rs115541804" "rs12044651" "rs11162246" ... $ MAF : num [1:4827] 0.4487 0.0288 0.4487 0.2715 0.0125 ... $ position: int [1:4827] 77445220 77445816 77446133 77446512 77446868 77447021 77447110 77447313 77447463 77447547 ... $ LD : num [1:4827, 1:4827] 1 -0.1833 0.996 -0.5917 -0.0828 ... ..- attr(*, "dimnames")=List of 2 .. ..$ : chr [1:4827] "rs12046188" "rs115541804" "rs12044651" "rs11162246" ... .. ..$ : chr [1:4827] "rs12046188" "rs115541804" "rs12044651" "rs11162246" ...
-
The code
S3= coloc::runsusie(coloc_data$gwas)
- The error
running max iterations: 100 WARNING: Providing the sample size (n), or even a rough estimate of n, is highly recommended. Without n, the implicit assumption is n is large (Inf) and the effect sizes are small (close to zero). WARNING: XtX is not symmetric; forcing XtX to be symmetric by replacing XtX with (XtX + t(XtX))/2 Error in susie_suff_stat(XtX = R, Xty = z, n = 2, yty = 1, scaled_prior_variance = prior_variance, : The estimated prior variance is unreasonably large. This is usually caused by mismatch between the summary statistics and the LD matrix. Please check the input. In addition: Warning message: In check_dataset(d, suffix, req = c("beta", "varbeta", "LD", "snp")) : minimum p value is: 4.0718e-06 If this is what you expected, this is not a problem. If this is not as small as you expected, please check the 02_data vignette.
- How to fix it
I think the problem in the code is this: https://github.com/chr1swallace/coloc/blob/27bed11a3e3fc22e258272b5929945770ab4b9ca/R/susie.R#L392https://eur03.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2Fchr1swallace%2Fcoloc%2Fblob%2F27bed11a3e3fc22e258272b5929945770ab4b9ca%2FR%2Fsusie.R%23L392&data=05%7C01%7Ccew54%40universityofcambridgecloud.onmicrosoft.com%7Ce1ff1fc6132b4a68e30c08dacbbfd2e6%7C49a50445bdfa4b79ade3547b4f3986e9%7C1%7C0%7C638046321980735828%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C&sdata=p%2FFclOUjnNmKJsXifHFABZqduXpGw2XHgaPdFnTAm0A%3D&reserved=0
susie_args=list(...) Error: '...' used in an incorrect context
Please, see an example how susie_rss works when I properly pass the sample size
4.1 Calling susie_rss
z and R are extracted from _colocdata$gwas same way you did in coloc::runsusie()
S3 = susieR::susie_rss(z = z, R = LD, n = n)
4.2. The output
WARNING: XtX is not symmetric; forcing XtX to be symmetric by replacing XtX with (XtX + t(XtX))/2 Warning message: In susie_suff_stat(XtX = XtX, Xty = Xty, n = n, yty = (n - 1) * : IBSS algorithm did not converge in 100 iterations! Please check consistency between summary statistics and LD matrix. See https://stephenslab.github.io/susieR/articles/susierss_diagnostic.html
str(S3) List of 18 $ alpha : num [1:10, 1:4827] 0 0 0 0 0 0 0 0 0 0 ... ..- attr(, "dimnames")=List of 2 .. ..$ : NULL .. ..$ : chr [1:4827] "rs12046188" "rs115541804" "rs12044651" "rs11162246" ... $ mu : num [1:10, 1:4827] -0.526 0.499 -0.531 0.501 -0.531 ... ..- attr(, "dimnames")=List of 2 .. ..$ : NULL .. ..$ : chr [1:4827] "rs12046188" "rs115541804" "rs12044651" "rs11162246" ... $ mu2 : num [1:10, 1:4827] 0.277 0.249 0.282 0.251 0.282 ... ..- attr(, "dimnames")=List of 2 .. ..$ : NULL .. ..$ : chr [1:4827] "rs12046188" "rs115541804" "rs12044651" "rs11162246" ... $ KL : num [1:10] 15.3 15.3 15.3 15.3 15.3 ... $ lbf : num [1:10] 379316 382738 386059 386440 386442 ... $ lbf_variable : num [1:10, 1:4827] 362 324 368 328 369 ... ..- attr(, "dimnames")=List of 2 .. ..$ : NULL .. ..$ : chr [1:4827] "rs12046188" "rs115541804" "rs12044651" "rs11162246" ... $ sigma2 : num 1 $ V : num [1:10] 285 287 290 290 290 ... $ pi : num [1:4827] 0.000207 0.000207 0.000207 0.000207 0.000207 ... $ null_index : num 0 $ XtXr : num [1:4827, 1] 6.784 -0.832 7.882 -4.908 -14.731 ... $ elbo : num [1:100] -3566 -3213 -2836 -2454 -2071 ... $ niter : int 100 $ converged : logi FALSE $ X_column_scale_factors: Named num [1:4827] 1 1 1 1 1 ... ..- attr(, "names")= chr [1:4827] "rs12046188" "rs115541804" "rs12044651" "rs11162246" ... $ intercept : num NA $ sets :List of 5 ..$ cs :List of 2 .. ..$ L1: int 1896 .. ..$ L2: int 1839 ..$ purity :'data.frame': 2 obs. of 3 variables: .. ..$ min.abs.corr : num [1:2] 1 1 .. ..$ mean.abs.corr : num [1:2] 1 1 .. ..$ median.abs.corr: num [1:2] 1 1 ..$ cs_index : int [1:2] 1 2 ..$ coverage : num [1:2] 1 1 ..$ requested_coverage: num 0.95 $ pip : Named num [1:4827] 0 0 0 0 0 0 0 0 0 0 ... ..- attr(, "names")= chr [1:4827] "rs12046188" "rs115541804" "rs12044651" "rs11162246" ...
- attr(*, "class")= chr "susie"
- sessionInfo()
sessionInfo() R version 4.2.0 (2022-04-22) Platform: x86_64-pc-linux-gnu (64-bit) Running under: Ubuntu 20.04.5 LTS
Matrix products: default BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/liblapack.so.3
locale: [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8 LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8 [7] LC_PAPER=en_US.UTF-8 LC_NAME=C LC_ADDRESS=C LC_TELEPHONE=C LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
attached base packages: [1] stats4 stats graphics grDevices utils datasets methods base
other attached packages: [1] rlang_1.0.5 susieR_0.12.28 snpStats_1.48.0 Matrix_1.4-1 [5] survival_3.3-1 SNPlocs.Hsapiens.dbSNP144.GRCh37_0.99.20 BSgenome_1.66.1 rtracklayer_1.58.0 [9] Biostrings_2.66.0 XVector_0.38.0 GenomicRanges_1.50.1 GenomeInfoDb_1.34.3 [13] IRanges_2.32.0 S4Vectors_0.36.0 BiocGenerics_0.44.0 data.table_1.14.2 [17] forcats_0.5.2 stringr_1.4.1 dplyr_1.0.10 purrr_0.3.4 [21] readr_2.1.2 tidyr_1.2.0 tibble_3.1.8 ggplot2_3.3.6 [25] tidyverse_1.3.2 reticulate_1.26 coloc_5.1.0 echolocatoR_2.0.3
loaded via a namespace (and not attached): [1] rappdirs_0.3.3 GGally_2.1.2 R.methodsS3_1.8.2 echoLD_0.99.8 bit64_4.0.5 knitr_1.40 [7] irlba_2.3.5 DelayedArray_0.24.0 R.utils_2.12.0 rpart_4.1.16 KEGGREST_1.38.0 RCurl_1.98-1.8 [13] AnnotationFilter_1.22.0 generics_0.1.3 GenomicFeatures_1.50.2 RSQLite_2.2.16 proxy_0.4-27 chron_2.3-57 [19] bit_4.0.4 tzdb_0.3.0 xml2_1.3.3 lubridate_1.8.0 SummarizedExperiment_1.28.0 assertthat_0.2.1 [25] viridis_0.6.2 gargle_1.2.0 xfun_0.32 hms_1.1.2 fansi_1.0.3 restfulr_0.0.15 [31] progress_1.2.2 dbplyr_2.2.1 readxl_1.4.1 Rgraphviz_2.41.1 igraph_1.3.4 DBI_1.1.3 [37] htmlwidgets_1.5.4 reshape_0.8.9 downloadR_0.99.5 googledrive_2.0.0 ellipsis_0.3.2 ggnewscale_0.4.7 [43] backports_1.4.1 biomaRt_2.54.0 deldir_1.0-6 MatrixGenerics_1.10.0 vctrs_0.4.1 Biobase_2.58.0 [49] ensembldb_2.22.0 cachem_1.0.6 withr_2.5.0 checkmate_2.1.0 GenomicAlignments_1.34.0 prettyunits_1.1.1 [55] cluster_2.1.3 ape_5.6-2 dir.expiry_1.6.0 lazyeval_0.2.2 crayon_1.5.1 basilisk.utils_1.10.0 [61] crul_1.2.0 pkgconfig_2.0.3 nlme_3.1-159 ProtGenerics_1.30.0 XGR_1.1.8 nnet_7.3-17 [67] pals_1.7 lifecycle_1.0.1 filelock_1.0.2 httpcode_0.3.0 BiocFileCache_2.6.0 modelr_0.1.9 [73] echotabix_0.99.8 dichromat_2.0-0.1 cellranger_1.1.0 matrixStats_0.62.0 graph_1.76.0 osfr_0.2.8 [79] boot_1.3-28 reprex_2.0.2 base64enc_0.1-3 googlesheets4_1.0.1 png_0.1-7 viridisLite_0.4.1 [85] rjson_0.2.21 rootSolve_1.8.2.3 bitops_1.0-7 R.oo_1.25.0 ggnetwork_0.5.10 blob_1.2.3 [91] mixsqp_0.3-43 echoplot_0.99.6 dnet_1.1.7 jpeg_0.1-9 echodata_0.99.16 scales_1.2.1 [97] memoise_2.0.1 magrittr_2.0.3 plyr_1.8.7 hexbin_1.28.2 zlibbioc_1.44.0 compiler_4.2.0 [103] echoconda_0.99.8 BiocIO_1.8.0 RColorBrewer_1.1-3 catalogueR_1.0.0 Rsamtools_2.14.0 cli_3.3.0 [109] echoannot_0.99.10 patchwork_1.1.2 htmlTable_2.4.1 Formula_1.2-4 MASS_7.3-58.1 tidyselect_1.1.2 [115] stringi_1.7.8 yaml_2.3.5 supraHex_1.35.0 latticeExtra_0.6-30 ggrepel_0.9.1 grid_4.2.0 [121] VariantAnnotation_1.44.0 tools_4.2.0 lmom_2.9 parallel_4.2.0 rstudioapi_0.14 foreign_0.8-82 [127] piggyback_0.1.3 gridExtra_2.3 gld_2.6.5 RcppZiggurat_0.1.6 digest_0.6.29 BiocManager_1.30.18 [133] proto_1.0.0 Rcpp_1.0.9 broom_1.0.1 OrganismDbi_1.40.0 httr_1.4.4 AnnotationDbi_1.60.0 [139] RCircos_1.2.2 ggbio_1.46.0 biovizBase_1.46.0 colorspace_2.0-3 rvest_1.0.3 XML_3.99-0.10 [145] fs_1.5.2 splines_4.2.0 RBGL_1.74.0 expm_0.999-6 echofinemap_0.99.4 basilisk_1.9.12 [151] Exact_3.1 mapproj_1.2.8 jsonlite_1.8.0 Rfast_2.0.6 R6_2.5.1 Hmisc_4.7-1 [157] gsubfn_0.7 pillar_1.8.1 htmltools_0.5.3 glue_1.6.2 fastmap_1.1.0 DT_0.24 [163] BiocParallel_1.32.1 class_7.3-20 codetools_0.2-18 maps_3.4.0 mvtnorm_1.1-3 utf8_1.2.2 [169] lattice_0.20-45 sqldf_0.4-11 curl_4.3.2 DescTools_0.99.46 zip_2.2.0 openxlsx_4.2.5 [175] interp_1.1-3 munsell_0.5.0 e1071_1.7-11 GenomeInfoDbData_1.2.9 haven_2.5.1 reshape2_1.4.4 [181] gtable_0.3.1
— Reply to this email directly, view it on GitHubhttps://eur03.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2Fchr1swallace%2Fcoloc%2Fissues%2F107&data=05%7C01%7Ccew54%40universityofcambridgecloud.onmicrosoft.com%7Ce1ff1fc6132b4a68e30c08dacbbfd2e6%7C49a50445bdfa4b79ade3547b4f3986e9%7C1%7C0%7C638046321980735828%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C&sdata=xqmJiDiroDG3W4exQuxq9az8t9fEQhprzL2bFn5HTQU%3D&reserved=0, or unsubscribehttps://eur03.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2Fnotifications%2Funsubscribe-auth%2FAAQWR2BO3L6UDMYXPP34J6TWJNWOHANCNFSM6AAAAAASGTVTKM&data=05%7C01%7Ccew54%40universityofcambridgecloud.onmicrosoft.com%7Ce1ff1fc6132b4a68e30c08dacbbfd2e6%7C49a50445bdfa4b79ade3547b4f3986e9%7C1%7C0%7C638046321980735828%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C&sdata=cz2YwK8VlM0Ot7d%2F7zD8aBxXV13EchaJtAOSo3lOHAs%3D&reserved=0. You are receiving this because you are subscribed to this thread.Message ID: @.***>
zy8281263
anbai106
pcarbo
Hi,
I just looked at the wrapper to call susie_rss().
Even though you add n as an argument to susie_args list, , and this should be properly passed to the do.call function wrapping susie_rss, when the function runs the n argument is missing.
1. The data
2. The code
3. The error
4. How to fix it
I think the problem in the code is this: https://github.com/chr1swallace/coloc/blob/27bed11a3e3fc22e258272b5929945770ab4b9ca/R/susie.R#L392
Please, see an example how susie_rss works when I properly pass the sample size
4.1 Calling susie_rss
4.2. The output
4. sessionInfo()