N argument missing when calling susieR::susie_rss()

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

 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" ...

2. The code

S3= coloc::runsusie(coloc_data$gwas)

3. 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.

4. How to fix it

I think the problem in the code is this:

coloc/R/susie.R

Line 392 in 27bed11

susie_args=list(...)

> 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 _coloc_data$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"

4. 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

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. 1. 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" ... 2. The code S3= coloc::runsusie(coloc_data$gwas) 3. 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. 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<https://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 _coloc_data$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" 4. 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 GitHub<https://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 unsubscribe<https://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: ***@***.***>

Hi,
I got the same issue with AMCalejandro due to N argument missing. The latest coloc is still 5.1.0.1 from CRAN, the issue was solved when downloading the latest coloc 5.2.1 from your github. Thanks!

thanks. I'll try to push a new version to CRAN.

…

-- https://chr1swallace.github.io<https://chr1swallace.github.io/>

________________________________ From: zy8281263 ***@***.***> Sent: Friday, May 12, 2023 8:38 AM To: chr1swallace/coloc ***@***.***> Cc: Chris Wallace ***@***.***>; Comment ***@***.***> Subject: Re: [chr1swallace/coloc] N argument missing when calling susieR::susie_rss() (Issue #107) Hi, I got the same issue with AMCalejandro due to N argument missing. The latest coloc is still 5.1.0.1 from CRAN, the issue was solved when downloading the latest coloc 5.2.1 from your github. Thanks! — Reply to this email directly, view it on GitHub<#107 (comment)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/AAQWR2DGKFSCLDEQ2WEXTJLXFXSGNANCNFSM6AAAAAASGTVTKM>. You are receiving this because you commented.Message ID: ***@***.***>

您好，我已收到您的邮件。我将尽快给您回复。

Hey,

I got a similar error while running "runsusie" from version 5.2.2:

5: 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

I don't see any problem with my LD matrix, which I used plink --bfile rs2744575_locus --r square --write-snplist to compute. In addition, the coloc.abf() gave reasonable results.

Any ideas?

Thanks

您好，我已收到您的邮件。我将尽快给您回复。

Can you check whether it converges in more iterations?

…

-- https://chr1swallace.github.io

________________________________ From: Junhao WEN ***@***.***> Sent: Thursday, June 8, 2023 9:21:40 PM To: chr1swallace/coloc ***@***.***> Cc: Chris Wallace ***@***.***>; Comment ***@***.***> Subject: Re: [chr1swallace/coloc] N argument missing when calling susieR::susie_rss() (Issue #107) Hey, I got a similar error while running "runsusie" from version 5.2.2: 5: 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 I don't see any problem with my LD matrix, which I used plink --bfile rs2744575_locus --r square --write-snplist to compute. In addition, the coloc.abf() gave reasonable results. Any ideas? Thanks — Reply to this email directly, view it on GitHub<#107 (comment)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/AAQWR2HXUPZP63QPTDZ6BVDXKIX5JANCNFSM6AAAAAASGTVTKM>. You are receiving this because you commented.Message ID: ***@***.***>

Same when I increased the iteration to be 1000 or 10000 Téléchargez Outlook pour iOS<https://aka.ms/o0ukef>

…

________________________________ De : Chris Wallace ***@***.***> Envoyé : Thursday, June 8, 2023 5:49:43 PM À : chr1swallace/coloc ***@***.***> Cc : Junhao WEN ***@***.***>; Comment ***@***.***> Objet : Re: [chr1swallace/coloc] N argument missing when calling susieR::susie_rss() (Issue #107) Can you check whether it converges in more iterations?

-- https://chr1swallace.github.io

________________________________ From: Junhao WEN ***@***.***> Sent: Thursday, June 8, 2023 9:21:40 PM To: chr1swallace/coloc ***@***.***> Cc: Chris Wallace ***@***.***>; Comment ***@***.***> Subject: Re: [chr1swallace/coloc] N argument missing when calling susieR::susie_rss() (Issue #107) Hey, I got a similar error while running "runsusie" from version 5.2.2: 5: 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 I don't see any problem with my LD matrix, which I used plink --bfile rs2744575_locus --r square --write-snplist to compute. In addition, the coloc.abf() gave reasonable results. Any ideas? Thanks — Reply to this email directly, view it on GitHub<#107 (comment)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/AAQWR2HXUPZP63QPTDZ6BVDXKIX5JANCNFSM6AAAAAASGTVTKM>. You are receiving this because you commented.Message ID: ***@***.***> — Reply to this email directly, view it on GitHub<#107 (comment)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/AD7X4KUDVAEAKBSNI2323KLXKJCHPANCNFSM6AAAAAASGTVTKM>. You are receiving this because you commented.Message ID: ***@***.***>

Did you check the consistency of the LD matrix as suggested? This is a Susie error, that can arise when the LD is calculated from a population that is not a close enough match to the GWAS sample or when alleles are not properly aligned.

…

-- https://chr1swallace.github.io

________________________________ From: Junhao WEN ***@***.***> Sent: Thursday, June 8, 2023 11:09:04 PM To: chr1swallace/coloc ***@***.***> Cc: Chris Wallace ***@***.***>; Comment ***@***.***> Subject: Re: [chr1swallace/coloc] N argument missing when calling susieR::susie_rss() (Issue #107) Same when I increased the iteration to be 1000 or 10000 Téléchargez Outlook pour iOS<https://aka.ms/o0ukef>

________________________________ De : Chris Wallace ***@***.***> Envoyé : Thursday, June 8, 2023 5:49:43 PM À : chr1swallace/coloc ***@***.***> Cc : Junhao WEN ***@***.***>; Comment ***@***.***> Objet : Re: [chr1swallace/coloc] N argument missing when calling susieR::susie_rss() (Issue #107) Can you check whether it converges in more iterations?

-- https://chr1swallace.github.io

________________________________ From: Junhao WEN ***@***.***> Sent: Thursday, June 8, 2023 9:21:40 PM To: chr1swallace/coloc ***@***.***> Cc: Chris Wallace ***@***.***>; Comment ***@***.***> Subject: Re: [chr1swallace/coloc] N argument missing when calling susieR::susie_rss() (Issue #107) Hey, I got a similar error while running "runsusie" from version 5.2.2: 5: 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 I don't see any problem with my LD matrix, which I used plink --bfile rs2744575_locus --r square --write-snplist to compute. In addition, the coloc.abf() gave reasonable results. Any ideas? Thanks — Reply to this email directly, view it on GitHub<#107 (comment)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/AAQWR2HXUPZP63QPTDZ6BVDXKIX5JANCNFSM6AAAAAASGTVTKM>. You are receiving this because you commented.Message ID: ***@***.***> — Reply to this email directly, view it on GitHub<#107 (comment)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/AD7X4KUDVAEAKBSNI2323KLXKJCHPANCNFSM6AAAAAASGTVTKM>. You are receiving this because you commented.Message ID: ***@***.***> — Reply to this email directly, view it on GitHub<#107 (comment)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/AAQWR2F5FAEGJJPDUQ6MF6LXKJEQBANCNFSM6AAAAAASGTVTKM>. You are receiving this because you commented.Message ID: ***@***.***>

I am checking now by following that page.

I doubt that is the reason, because I am using local genotype data to run GWAS, not GWAS summary stats from previous literature.

I will be back -:

I got some diagnostic results by following that post (https://stephenslab.github.io/susieR/articles/susierss_diagnostic.html) using my own data (Section: LD information from the original genotype data).

Here is the regional Manhattan plot for this locus (from the coloc.abf function):

Here is the plot for the z-scores (Y-axis is the z-score, instead of the -log10(P) in the figure above):

The estimated s explain the consistency between the z scores, and the LD matrix is 0.02832 (s = estimate_s_rss(z_scores, Rin, n=111301)). Their simulation data gave a much smaller s value, compared to mine.

Note: the LD matrix (Rin, calculated by Plink using the --r argument) gave this warning:

WARNING: The matrix R is not positive semidefinite. Negative eigenvalues are set to zero.

Finally, the distribution of the z-scores is here:

I am new to fine mapping and colocalization. Can you please help to give some insights on these?
@pcarbo

Thanks for this wonderful tool!

@anbai106 Overall, these results look good. (You can ignore the warning—not important.) Potentially there are a few SNPs you could remove judging by the z-score scatterplot, but doesn't seem like a major concern. It looks like you are okay to move to the next step and run susie_rss by following the susieR vignettes, then examing the susie_rss outputs closely to see if they make sense. (Or, if your aim is to perform colocalization, you may want to run susie_rss separately first and check that the results make sense before runnning coloc with susie.)

@pcarbo Thank you for your reply!

susieR version (susieR_0.12.35) on Ubuntu 18.04

So, now I am following the vignettes from susieR to do fine-mapping with summary statistics (https://stephenslab.github.io/susieR/articles/finemapping_summary_statistics.html):

My code is here:

library(susieR)
library(curl)
rm(list = ls()) ## remove global environment

musculo_tsv <- "rs2744575_coloc_prepare.tsv"
musculo_ld <- "plink.ld" #(in-sample LD via plink --r square)
musculo_ld_snp <- "plink.snplist"
musculo_type <- "rs2744575_pheno_type.txt"
musculo_pheno_std <- "rs2744575_pheno_std.txt"
df_data_musculo <- read.table(musculo_tsv, header=T, sep='\t', quote="")
df_data_musculo_ld <- read.table(musculo_ld, header=F, sep='\t', quote="")
df_data_musculo_ld_snp <- read.table(musculo_ld_snp, header=F, sep='\t', quote="")
df_data_musculo_type <- read.table(musculo_type, header=F, sep='\t', quote="")
df_data_musculo_std <- read.table(musculo_pheno_std, header=F, sep='\t', quote="")
D_musculo = as.list(df_data_musculo)
D_musculo$sebeta <- D_musculo$varbeta
D_musculo$varbeta = D_musculo$varbeta^2
D_musculo$sdY = df_data_musculo_std$V1
D_musculo$N = 111301
D_musculo$type = df_data_musculo_type$V1
colnames(df_data_musculo_ld) <- df_data_musculo_ld_snp$V1
rownames(df_data_musculo_ld) <- df_data_musculo_ld_snp$V1
D_musculo$LD = data.matrix(df_data_musculo_ld)
check_dataset(D_musculo)

z_scores <- D_musculo$beta / D_musculo$sebeta
Rin = D_musculo$LD
attr(Rin, "eigen") = eigen(Rin, symmetric = TRUE)
susie_plot(z_scores, y = "z")
s = estimate_s_rss(z_scores, Rin, n=111301)
condz_in = kriging_rss(z_scores, Rin, n=111301)
condz_in$plot

fitted_rss1 <- susie_rss(bhat=D_musculo$beta, shat=D_musculo$sebeta, n=D_musculo$N , R=D_musculo$LD, var_y=D_musculo$sdY, L=10,
estimate_residual_variance = TRUE)
summary(fitted_rss1)$cs

This gave the following results: There are many cs with cs_log10bf = Inf.

cs cs_log10bf cs_avg_r2 cs_min_r2 variable
1 5 Inf 1 1 1193,1195
2 1 212.8518 1 1 1497
3 2 Inf 1 1 1256
4 3 Inf 1 1 1185
5 4 Inf 1 1 1373
6 6 Inf 1 1 1265
7 7 Inf 1 1 1374
8 8 Inf 1 1 790
9 9 Inf 1 1 1120
10 10 Inf 1 1 1180,1181

So I found some other people suggest that this may be due to the convergence issues:

fitted_rss1 <- susie_rss(bhat=D_musculo$beta, shat=D_musculo$sebeta, n=D_musculo$N , R=D_musculo$LD, var_y=D_musculo$sdY, L=10,
estimate_residual_variance = TRUE, max_iter = 250, verbose = TRUE)

Then I got errors:

[1] "objective:-157779.641408244"
[1] "objective:-157779.587498904"
[1] "objective:-157762.674670233"
[1] "objective:-157762.671285953"
[1] "objective:-157756.556927372"
[1] "objective:-157756.556291964"
[1] "objective:-157753.422088498"
[1] "objective:-157753.421882784"
[1] "objective:-157750.32962681"
[1] "objective:-157750.329439944"
[1] "objective:-157745.717673608"
[1] "objective:-157745.717257631"
[1] "objective:-157739.573213663"
[1] "objective:-157739.572570525"
[1] "objective:-157733.589148574"
[1] "objective:-157733.588670104"
[1] "objective:-157727.09907281"
[1] "objective:-157727.098497427"
[1] "objective:-157720.73753681"
[1] "objective:-157720.736921909"
[1] "objective:-157715.002562835"
[1] "objective:-157715.002074037"
[1] "objective:-157709.462236342"
[1] "objective:-157709.46180983"
[1] "objective:-157703.857401504"
[1] "objective:-157703.856961682"
[1] "objective:-157697.553902515"
[1] "objective:-157697.553289637"
[1] "objective:-157691.141516755"
[1] "objective:-157691.140954493"
[1] "objective:-157684.626734458"
[1] "objective:-157684.626182029"
[1] "objective:-157677.715920462"
[1] "objective:-157677.715340372"
[1] "objective:-157670.456315951"
[1] "objective:-157670.455717035"
[1] "objective:-157662.989589638"
[1] "objective:-157662.988995681"
[1] "objective:-157655.406484043"
[1] "objective:-157655.405894388"
[1] "objective:-157647.752420657"
[1] "objective:-157647.751830875"
[1] "objective:-157640.050415887"
[1] "objective:-157640.049823363"
[1] "objective:-157632.313645501"
[1] "objective:-157632.313047182"
[1] "objective:-157624.546299243"
[1] "objective:-157624.545685175"
[1] "objective:-157616.719765236"
[1] "objective:-157616.719088297"
[1] "objective:-157608.433981795"
[1] "objective:-157608.432761337"
[1] "objective:-157595.768309387"
[1] "objective:-157595.765177119"
[1] "objective:-157583.168103205"
[1] "objective:-157583.166448337"
[1] "objective:-157571.124442387"
[1] "objective:-157571.123052082"
[1] "objective:-157559.759963033"
[1] "objective:-157559.758623324"
[1] "objective:-157548.632169724"
[1] "objective:-157548.630962609"
[1] "objective:-157538.281497756"
[1] "objective:-157538.280901973"
[1] "objective:-157522.713168083"
[1] "objective:-157522.709388463"
[1] "objective:-157507.901306927"
[1] "objective:-157507.897984015"
[1] "objective:-157492.479520809"
[1] "objective:-157492.476322775"
[1] "objective:-157477.04362092"
[1] "objective:-157477.040939776"
[1] "objective:-157461.475651436"
[1] "objective:-157461.472864212"
[1] "objective:-157445.092150976"
[1] "objective:-157445.08814025"
[1] "objective:-157429.844624384"
[1] "objective:-157429.842177463"
[1] "objective:-157413.969939879"
[1] "objective:-157413.967536107"
[1] "objective:-157397.736816515"
[1] "objective:-157397.734296245"
[1] "objective:-157380.854705906"
[1] "objective:-157380.851808495"
[1] "objective:-157357.992359466"
[1] "objective:-157357.98761905"
[1] "objective:-157330.55255297"
[1] "objective:-157330.545328654"
[1] "objective:-157295.065213535"
[1] "objective:-157295.053047382"
[1] "objective:-157260.201722552"
[1] "objective:-157260.189955514"
[1] "objective:-157222.549227981"
[1] "objective:-157222.535880068"
[1] "objective:-157182.805992707"
[1] "objective:-157182.791095907"
[1] "objective:-157140.927480801"
[1] "objective:-157140.910991582"
[1] "objective:-157096.837805707"
[1] "objective:-157096.8196853"
[1] "objective:-157050.488219843"
[1] "objective:-157050.468341442"
[1] "objective:-157001.798542246"
[1] "objective:-157001.776709014"
[1] "objective:-156950.672046548"
[1] "objective:-156950.648046014"
[1] "objective:-156897.00602622"
[1] "objective:-156896.9796414"
[1] "objective:-156840.693983624"
[1] "objective:-156840.664992549"
[1] "objective:-156781.625239899"
[1] "objective:-156781.593419287"
[1] "objective:-156719.704686298"
[1] "objective:-156719.669864446"
[1] "objective:-156655.191211361"
[1] "objective:-156655.153695154"
[1] "objective:-156581.629483863"
[1] "objective:-156581.579396382"
[1] "objective:-156495.10138558"
[1] "objective:-156495.030892854"
[1] "objective:-156413.440621029"
[1] "objective:-156413.37803353"
[1] "objective:-156328.806118123"
[1] "objective:-156328.740624731"
[1] "objective:-156240.399389765"
[1] "objective:-156240.328401606"
[1] "objective:-156147.844267634"
[1] "objective:-156147.766592481"
[1] "objective:-156050.881664304"
[1] "objective:-156050.796482814"
[1] "objective:-155949.250016236"
[1] "objective:-155949.156484338"
[1] "objective:-155842.681418566"
[1] "objective:-155842.578626426"
[1] "objective:-155730.901946163"
[1] "objective:-155730.788905304"
[1] "objective:-155613.630011611"
[1] "objective:-155613.505642693"
[1] "objective:-155490.573979964"
[1] "objective:-155490.437099376"
[1] "objective:-155361.429805513"
[1] "objective:-155361.27911097"
[1] "objective:-155225.879123303"
[1] "objective:-155225.713179106"
[1] "objective:-155083.592275371"
[1] "objective:-155083.409510034"
[1] "objective:-154934.316403697"
[1] "objective:-154934.115353057"
[1] "objective:-154780.012428335"
[1] "objective:-154779.797957557"
[1] "objective:-154625.378501939"
[1] "objective:-154625.160785543"
[1] "objective:-154449.466089356"
[1] "objective:-154449.186172508"
[1] "objective:-154264.379014184"
[1] "objective:-154264.070069616"
[1] "objective:-154069.289505886"
[1] "objective:-154068.946529728"
[1] "objective:-153863.736795472"
[1] "objective:-153863.356214806"
[1] "objective:-153647.22552428"
[1] "objective:-153646.803433587"
[1] "objective:-153419.199821073"
[1] "objective:-153418.731797877"
[1] "objective:-153179.039370098"
[1] "objective:-153178.520391818"
[1] "objective:-152926.047980092"
[1] "objective:-152925.472311162"
[1] "objective:-152659.416493096"
[1] "objective:-152658.777447493"
[1] "objective:-152378.107903133"
[1] "objective:-152377.397168716"
[1] "objective:-152080.511521655"
[1] "objective:-152079.716869553"
[1] "objective:-151763.995585352"
[1] "objective:-151763.095945484"
[1] "objective:-151427.438644881"
[1] "objective:-151426.41893644"
[1] "objective:-151072.018343978"
[1] "objective:-151070.883343971"
[1] "objective:-150697.283257534"
[1] "objective:-150696.021662662"
[1] "objective:-150301.623883441"
[1] "objective:-150300.217470743"
[1] "objective:-149883.472534311"
[1] "objective:-149881.901453721"
[1] "objective:-149441.281202256"
[1] "objective:-149439.523799657"
[1] "objective:-148973.434582822"
[1] "objective:-148971.467098034"
[1] "objective:-148478.187066108"
[1] "objective:-148475.982806644"
[1] "objective:-147953.635352173"
[1] "objective:-147951.163498303"
[1] "objective:-147397.723864843"
[1] "objective:-147394.948745279"
[1] "objective:-146808.23876557"
[1] "objective:-146805.119375956"
[1] "objective:-146182.782223391"
[1] "objective:-146179.271483824"
[1] "objective:-145518.740319609"
[1] "objective:-145514.784005242"
[1] "objective:-144813.250623868"
[1] "objective:-144808.786027915"
[1] "objective:-144063.168724972"
[1] "objective:-144058.123053157"
[1] "objective:-143265.031422787"
[1] "objective:-143259.319863475"
[1] "objective:-142415.014472697"
[1] "objective:-142408.537853832"
[1] "objective:-141508.883088179"
[1] "objective:-141501.525005225"
[1] "objective:-140541.933621155"
[1] "objective:-140533.556891829"
[1] "objective:-139508.925954246"
[1] "objective:-139499.368240482"
[1] "objective:-138404.016876414"
[1] "objective:-138393.085497924"
[1] "objective:-137220.832764042"
[1] "objective:-137208.301585651"
[1] "objective:-135957.058595322"
[1] "objective:-135942.769984841"
[1] "objective:-134539.614006501"
[1] "objective:-134521.609413301"
[1] "objective:-132502.986575204"
[1] "objective:-132465.431546825"
[1] "objective:-130472.339805861"
[1] "objective:-130435.745798919"
[1] "objective:-128385.77627683"
[1] "objective:-128347.040957983"
[1] "objective:-126184.349046833"
[1] "objective:-126141.179712989"
[1] "objective:-123820.9205687"
[1] "objective:-123771.137321276"
[1] "objective:-121251.580312435"
[1] "objective:-121192.734570339"
[1] "objective:-118430.576304442"
[1] "objective:-118359.641136479"
[1] "objective:-115306.381951108"
[1] "objective:-115219.393566969"
[1] "objective:-111817.686393426"
[1] "objective:-111709.238648589"
[1] "objective:-107888.406661904"
[1] "objective:-107750.864078421"
[1] "objective:-103420.672448737"
[1] "objective:-103242.865539494"
[1] "objective:-98284.1418798462"
[1] "objective:-98049.0806464103"
[1] "objective:-92298.6602874564"
[1] "objective:-91979.2767487976"
[1] "objective:-85204.2728749045"
[1] "objective:-84754.8951237346"
[1] "objective:-76605.4453367068"
[1] "objective:-75943.0739354484"
[1] "objective:-65857.3731577971"
[1] "objective:-64815.2037660855"
[1] "objective:-51804.1296469366"
[1] "objective:-49994.4260833244"
[1] "objective:-32060.0257630374"
[1] "objective:-28345.384686537"
[1] "objective:-401.049625791777"
[1] "objective:10469.5417334626"
[1] "objective:69736.492225531"
Error in susie_suff_stat(XtX = XtX, Xty = Xty, n = n, yty = (n - 1) * :
Estimating residual variance failed: the estimated value is negative

Let me know if you want me to share the data.

@chr1swallace, so it seems that the data have no issues. Why does the runsusie function give the errors? How can I proceed if I want to use this function, instead of using the original susie package?

runsusie is just a wrapper around susie_rss that 1. accepts data in the same format as other coloc functions 2. calls susie_rss 3. annotates the output of susie_rss with snp identifiers so that the results can be aligned with results from another trait The error you are getting is from susie_rss, so if you can really call that function directly without issues, I would like to see your data to understand what runsusie is doing to cause the issue. For now, I have separated step 3 above, so that you could call susie_rss directly, then run annotate_susie, and pass the result to coloc_susie. ie res = susie_rss(...) res = annotate_susie(res, snpids) coloc_susie(res, [other dataset]) This is available in the latest github version of coloc.

…

-- https://chr1swallace.github.io<https://chr1swallace.github.io/>

________________________________ From: Junhao WEN ***@***.***> Sent: Wednesday, June 14, 2023 1:57 AM To: chr1swallace/coloc ***@***.***> Cc: Chris Wallace ***@***.***>; Mention ***@***.***> Subject: Re: [chr1swallace/coloc] N argument missing when calling susieR::susie_rss() (Issue #107) @chr1swallace<https://github.com/chr1swallace>, so it seems that the data have no issues. Why does the runsusie function give the errors? How can I proceed if I want to use this function, instead of using the original susie package? — Reply to this email directly, view it on GitHub<#107 (comment)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/AAQWR2B5V4BVFOACCROQSODXLEEAHANCNFSM6AAAAAASGTVTKM>. You are receiving this because you were mentioned.Message ID: ***@***.***>

@anbai106 Regarding the susie error, try again by settiing estimate_residual_variance = FALSE.

We have also found that setting estimate_prior_method = "EM" sometimes helps (somewhat surprisingly—we don't quite understand why).

@pcarbo

I saw the suggestion that you made here: https://stephenslab.github.io/susieR/articles/finemapping_summary_statistics.html

To confirm that I am using the in-sample LD, which in the documentation suggests setting estimate_residual_variance = TRUE.

So I reran the function by following your suggestion 1 and 2:

Suggestion 1:

fitted_rss1 <- susie_rss(bhat=D_musculo$beta, shat=D_musculo$sebeta, n=D_musculo$N , R=D_musculo$LD, var_y=D_musculo$sdY, L=10,
                         estimate_residual_variance = FALSE, verbose = TRUE)
summary(fitted_rss1)$cs

The output is:

[1] "objective:-157779.641408244"
[1] "objective:-157762.756068947"
[1] "objective:-157756.657349565"
[1] "objective:-157753.534467443"
[1] "objective:-157750.465269108"
[1] "objective:-157745.883921601"
[1] "objective:-157739.760616217"
[1] "objective:-157733.787467113"
[1] "objective:-157727.309659819"
[1] "objective:-157720.976469878"
[1] "objective:-157715.263702416"
[1] "objective:-157709.746406992"
[1] "objective:-157704.169239631"
[1] "objective:-157697.901244172"
[1] "objective:-157691.522897711"
[1] "objective:-157685.034423654"
[1] "objective:-157678.156477157"
[1] "objective:-157670.931608979"
[1] "objective:-157663.500013296"
[1] "objective:-157655.952550686"
[1] "objective:-157648.335027928"
[1] "objective:-157640.670702673"
[1] "objective:-157632.972964069"
[1] "objective:-157625.24662591"
[1] "objective:-157617.466883113"
[1] "objective:-157609.299841714"
[1] "objective:-157596.966121912"
[1] "objective:-157584.378838127"
[1] "objective:-157572.360937088"
[1] "objective:-157561.008666021"
[1] "objective:-157549.956957514"
[1] "objective:-157539.447542277"
[1] "objective:-157524.391320436"
[1] "objective:-157509.759960009"
[1] "objective:-157494.498971446"
[1] "objective:-157479.197826273"
[1] "objective:-157463.793321649"
[1] "objective:-157447.651979678"
[1] "objective:-157432.371230713"
[1] "objective:-157416.731189028"
[1] "objective:-157400.658162667"
[1] "objective:-157384.206354768"
[1] "objective:-157362.611413473"
[1] "objective:-157336.725499403"
[1] "objective:-157302.394860125"
[1] "objective:-157267.913350862"
[1] "objective:-157230.897317032"
[1] "objective:-157191.876221454"
[1] "objective:-157150.807780167"
[1] "objective:-157107.601572863"
[1] "objective:-157062.213282006"
[1] "objective:-157014.572128645"
[1] "objective:-156964.590223314"
[1] "objective:-156912.173520299"
[1] "objective:-156857.224578855"
[1] "objective:-156799.642347501"
[1] "objective:-156739.339036044"
[1] "objective:-156676.532551273"
[1] "objective:-156605.71448285"
[1] "objective:-156521.490155557"
[1] "objective:-156442.414347904"
[1] "objective:-156360.555766032"
[1] "objective:-156275.186176614"
[1] "objective:-156185.956695626"
[1] "objective:-156092.637671841"
[1] "objective:-155995.000022114"
[1] "objective:-155892.811364008"
[1] "objective:-155785.836577073"
[1] "objective:-155673.83652284"
[1] "objective:-155556.566075443"
[1] "objective:-155433.772219699"
[1] "objective:-155305.192747162"
[1] "objective:-155170.560991697"
[1] "objective:-155029.703277068"
[1] "objective:-154884.659396797"
[1] "objective:-154739.481784485"
[1] "objective:-154574.803541293"
[1] "objective:-154402.160927365"
[1] "objective:-154220.81230122"
[1] "objective:-154030.422812368"
[1] "objective:-153830.637450693"
[1] "objective:-153621.058065532"
[1] "objective:-153401.239779101"
[1] "objective:-153170.680522386"
[1] "objective:-152928.787147765"
[1] "objective:-152674.776482248"
[1] "objective:-152407.414696221"
[1] "objective:-152124.805170005"
[1] "objective:-151826.244343593"
[1] "objective:-151512.816683187"
[1] "objective:-151184.447850685"
[1] "objective:-150840.129973292"
[1] "objective:-150478.890455208"
[1] "objective:-150099.820123892"
[1] "objective:-149702.010146651"
[1] "objective:-149284.517992594"
[1] "objective:-148846.351925957"
[1] "objective:-148386.474761265"
[1] "objective:-147903.80970306"
[1] "objective:-147397.237915452"
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
Browse[2]> summary(fitted_rss1)$cs
cs cs_log10bf cs_avg_r2 cs_min_r2 variable
1 5 Inf 1 1 1193,1195
2 1 172.2678 1 1 1497
3 2 Inf 1 1 1256
4 3 Inf 1 1 1185
5 4 Inf 1 1 1373
6 6 Inf 1 1 1265
7 7 Inf 1 1 1374
8 8 Inf 1 1 790
9 9 Inf 1 1 1120
10 10 Inf 1 1 1180,1181

Suggestion 2:

fitted_rss1 <- susie_rss(bhat=D_musculo$beta, shat=D_musculo$sebeta, n=D_musculo$N , R=D_musculo$LD, var_y=D_musculo$sdY, L=10,
                         estimate_residual_variance = FALSE, estimate_prior_method = "EM", verbose = TRUE)
summary(fitted_rss1)$cs

The output is:

[1] "objective:-157814.881509146"
[1] "objective:-157765.518602674"
[1] "objective:-157759.989527055"
[1] "objective:-157755.456909742"
[1] "objective:-157752.695451689"
[1] "objective:-157749.644145321"
[1] "objective:-157744.535190621"
[1] "objective:-157738.554966968"
[1] "objective:-157732.190853185"
[1] "objective:-157725.837635829"
[1] "objective:-157720.087468464"
[1] "objective:-157714.81439392"
[1] "objective:-157709.922462254"
[1] "objective:-157705.132141293"
[1] "objective:-157699.982567725"
[1] "objective:-157694.183839972"
[1] "objective:-157687.895900036"
[1] "objective:-157681.191929984"
[1] "objective:-157674.149855358"
[1] "objective:-157666.875192079"
[1] "objective:-157659.455692131"
[1] "objective:-157651.945275281"
[1] "objective:-157644.374283762"
[1] "objective:-157636.762940935"
[1] "objective:-157629.127990793"
[1] "objective:-157621.484455996"
[1] "objective:-157613.845664414"
[1] "objective:-157606.223058642"
[1] "objective:-157598.626189009"
[1] "objective:-157591.062875095"
[1] "objective:-157583.539442665"
[1] "objective:-157576.060967589"
[1] "objective:-157568.631493162"
[1] "objective:-157561.254210338"
[1] "objective:-157553.931601788"
[1] "objective:-157546.665555071"
[1] "objective:-157539.457450871"
[1] "objective:-157532.308231619"
[1] "objective:-157525.218454788"
[1] "objective:-157518.188334258"
[1] "objective:-157511.217772507"
[1] "objective:-157504.306385955"
[1] "objective:-157497.453525591"
[1] "objective:-157490.658294978"
[1] "objective:-157483.91956781"
[1] "objective:-157477.236007332"
[1] "objective:-157470.606089994"
[1] "objective:-157464.028135649"
[1] "objective:-157457.500346115"
[1] "objective:-157451.020853132"
[1] "objective:-157444.587775303"
[1] "objective:-157438.199281871"
[1] "objective:-157431.853659176"
[1] "objective:-157425.549374052"
[1] "objective:-157419.285127518"
[1] "objective:-157413.05989261"
[1] "objective:-157406.872931923"
[1] "objective:-157400.723793291"
[1] "objective:-157394.612285213"
[1] "objective:-157388.538436396"
[1] "objective:-157382.502445448"
[1] "objective:-157376.504627161"
[1] "objective:-157370.545360973"
[1] "objective:-157364.625045731"
[1] "objective:-157358.744063093"
[1] "objective:-157352.90275033"
[1] "objective:-157347.101382129"
[1] "objective:-157341.340160241"
[1] "objective:-157335.619209488"
[1] "objective:-157329.938578598"
[1] "objective:-157324.298244467"
[1] "objective:-157318.69811866"
[1] "objective:-157313.138055197"
[1] "objective:-157307.617858923"
[1] "objective:-157302.137293916"
[1] "objective:-157296.696091589"
[1] "objective:-157291.293958242"
[1] "objective:-157285.93058191"
[1] "objective:-157280.605638442"
[1] "objective:-157275.318796778"
[1] "objective:-157270.069723424"
[1] "objective:-157264.858086167"
[1] "objective:-157259.683557062"
[1] "objective:-157254.545814764"
[1] "objective:-157249.444546243"
[1] "objective:-157244.379447979"
[1] "objective:-157239.350226677"
[1] "objective:-157234.356599586"
[1] "objective:-157229.398294483"
[1] "objective:-157224.475049386"
[1] "objective:-157219.586612051"
[1] "objective:-157214.732739316"
[1] "objective:-157209.913196333"
[1] "objective:-157205.127755734"
[1] "objective:-157200.376196767"
[1] "objective:-157195.658304443"
[1] "objective:-157190.973868698"
[1] "objective:-157186.322683607"
[1] "objective:-157181.704546666"
[1] "objective:-157177.119258134"
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
Browse[2]> summary(fitted_rss1)$cs
cs cs_log10bf cs_avg_r2 cs_min_r2 variable
1 5 Inf 1.0000000 1.000000 1193,1195
2 2 Inf 1.0000000 1.000000 1256
3 3 Inf 1.0000000 1.000000 1185
4 4 Inf 0.9998408 0.999602 1192,1196,1206,1207,1223
5 1 34.96028 0.9070886 0.861364 1276,1296,1301

They gave quite different results for the non-INF cs_log10bf. The INF means that the model do not converge right? The second suggestion resulted in only 5 cs.

If I increased the max_iter = 500 for both suggestions 1 and 2:

Suggestion 1:

Error in susie_suff_stat(XtX = XtX, Xty = Xty, n = n, yty = (n - 1) * :
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.

Suggestion 2:

Warning message:
In susie_suff_stat(XtX = XtX, Xty = Xty, n = n, yty = (n - 1) * :
IBSS algorithm did not converge in 500 iterations!
Please check consistency between summary statistics and LD matrix.
See https://stephenslab.github.io/susieR/articles/susierss_diagnostic.html
Browse[2]> summary(fitted_rss1)$cs
cs cs_log10bf cs_avg_r2 cs_min_r2 variable
1 1 Inf 1 1 1374
2 2 Inf 1 1 1175
3 3 Inf 1 1 1185
4 4 Inf 1 1 1223
5 5 Inf 1 1 1195

Thanks for your help!

Best regards

@chr1swallace Thank you for your guidance, Chris.

I think I will run my analyses with coloc.abf for now. I got some reasonable results, but would like to ask some advice to interpret them.

The results summary is:

nsnps 1549
PP.H0.abf 2.27180411599461e-10
PP.H1.abf 0.117176654233345
PP.H2.abf 3.52382045115394e-11
PP.H3.abf 0.0173098867739821
PP.H4.abf 0.865513458730256

I found potential colocalization between my two traits under the H4 hypothesis threshold (>0.8). The sensitivity check is here:

rl)

My interpretation is as below (please correct me if I am wrong):

Bayesian colocalization analyses supported a potential causal association of SNPs within this locus with trait1 and trait2. The results showed a posterior possibility (PP) of one shared causal variant (PP.H4.ABF=0.8655) associated with both traits in the XXX mapped gene.

I have additional questions:

In the Posterior Probability plot, the green shadow regions indicate the colocalization results are not very robust - I must set the prior probability of p12 in the green region? I set this prior prob by default.
Since the dotted grey line (results) is inside the green region, the results are reasonable and acceptable?
As we can see, the locus on trait 2 does not achieve the genome-wide association threshold (-logP>7.3), does this imply the H1 hypothesis PP is 0.11? I read your paper, and you suggest that we perform colocalization only when both traits are significantly associated with both traits. Is the current situation OK? I actually choose the locus based on trait1 GWAS results.
How to choose the sensitivity and causal PP threshold? I saw some paper uses >0.9
For the colored SNP on the left panel for both trait1 (upper) and trait2 (lower), the darkest SNP has the highest probability of being the only causal SNP in the locus under the H4 hypothesis?

One more question, I also did fine-mapping using finemap.abf, how do we choose the threshold for SNP.PP to decide if the SNP is causal? Normally, the top lead SNP (most significant in the locus) has the highest SNP.PP value, but not often achieved the 0.8 thresholds.

Thanks for your help!

@anbai106 A log10 BF of Inf means that the support for the CS is very large. Separately, it does also look like the IBSS model fitting algorithm is having difficulty converging to a solution. I'm not sure why—it could be because it is a more challenging fine-mapping/colocalization problem. Or there could be some issues with your data.

@pcarbo
Thank you for your reply, I will double-check my own data.

But from the results that I got from coloc.abf, the results do make sense and correspond very well with my genetic correlation results between the two traits of interest. The only possibility of the problem of my data is the LD matrix...