chr1swallace/coloc

runsusie and susieR results differ

maguileraf opened this issue · 8 comments

maguileraf commented

Hi Chris,

I was using runsusie but for some regions, it doesn't give me any credible sets. However, when I run the same dataset using susieR, I get credible sets. Also, I compared the output from runsusie and susieR and they are different. Do you know why this is happening?

chr1swallace commented
maguileraf commented

Unfortunately, I can't share the data. I have attached an example of the output for both commands. Also, these are the commands I am using:
S1 <- runsusie(D1, estimate_residual_variance = F, n = 555634, refine = T)
gwas_susie <- susie_rss(gwas_for_coloc$Z, R = loci_r_matrix, L = 10, n = 555634,
estimate_residual_variance = F, refine=T)
They are both using the same LD matrix and same data. However, maybe the Z scores are being calculated differently? I have the latest version available in GitHub of susieR, maybe this could also make a difference.

Screenshot 2023-11-28 at 6 24 08 AM Screenshot 2023-11-28 at 6 24 26 AM
chr1swallace commented
maguileraf commented

I checked the z-scores and they are the same with my formula and yours. Attached is the distribution of the z-scores.
911736_990053zscore_hist_1112823.pdf

catheriz commented

This happened to me as well:

gwas_S=runsusie(gwas,n = unique(gwas$N),L = 10, max_iter = 2500,verbose=TRUE,coverage = 0.99, min_abs_corr = sqrt(0.4),refine = T)
running max iterations: 2500
[1] "objective:-16969.693590593"
[1] "objective:-16969.6211774467"
[1] "objective:-16969.5931406855"
[1] "objective:-16969.5784032943"
[1] "objective:-16969.5687705576"
[1] "objective:-16969.5610090115"
[1] "objective:-16969.5536898941"
[1] "objective:-16969.5459284366"
[1] "objective:-16969.5364229281"
[1] "objective:-16969.5228911058"
[1] "objective:-16969.5017541661"
[1] "objective:-16969.4676286733"
[1] "objective:-16969.4152138044"
[1] "objective:-16969.3516179333"
[1] "objective:-16969.3016357743"
[1] "objective:-16969.2752213206"
[1] "objective:-16969.2639810721"
[1] "objective:-16969.2594363808"
[1] "objective:-16969.2573272553"
[1] "objective:-16969.2560825934"
[1] "objective:-16969.2552228556"
converged: TRUE
summary(gwas_S)

Variables in credible sets:

variable variable_prob cs
354 1.0000000000 4
801 1.0000000000 3
838 0.1766969636 6
836 0.1736468743 6
810 0.1560065758 6
107 0.1044441996 5
104 0.1013487080 5
809 0.0925986824 6
81 0.0827069134 5
100 0.0581785893 5
852 0.0465301371 6
805 0.0302825176 6
844 0.0283636475 6
824 0.0281428096 6
848 0.0281107760 6
825 0.0280476100 6
826 0.0277783989 6
829 0.0262627277 6
818 0.0262582751 6
840 0.0257368057 6
813 0.0255036182 6
821 0.0223707009 6
874 0.0191023850 6
865 0.0188001994 6
138 0.0165209232 5
72 0.0163644076 5
58 0.0163203527 5
71 0.0161607759 5
68 0.0152525246 5
36 0.0151258413 5
24 0.0150888761 5
18 0.0150728397 5
15 0.0150464225 5
60 0.0149663343 5
76 0.0149004323 5
89 0.0148920532 5
75 0.0148654132 5
12 0.0148394278 5
184 0.0148378343 5
59 0.0148324238 5
35 0.0148087350 5
44 0.0148059259 5
64 0.0147951808 5
34 0.0147619432 5
43 0.0147575432 5
25 0.0147447166 5
28 0.0147048548 5
27 0.0147026981 5
39 0.0146483521 5
66 0.0146307209 5
54 0.0145901830 5
192 0.0145680505 5
30 0.0144262035 5
70 0.0143358401 5
29 0.0142752598 5
129 0.0140407037 5
172 0.0139046679 5
146 0.0138850559 5
47 0.0138343094 5
167 0.0138274744 5
152 0.0137749771 5
147 0.0137668958 5
38 0.0137462694 5
145 0.0136617665 5
63 0.0134298156 5
804 0.0131904877 6
141 0.0129086439 5
163 0.0118078825 5
96 0.0079051393 5
97 0.0048937345 5
491 0.0008941885 5
505 0.0008890025 5
473 0.0008244679 5
3 0.0004672167 5
497 0.0004122179 5

Credible sets summary:

cs cs_log10bf cs_avg_r2 cs_min_r2
3 16.227200 1.0000000 1.0000000
4 72.959516 1.0000000 1.0000000
5 44.486494 0.9780216 0.7662516
6 3.679031 0.8794713 0.6670446
variable
801
354
3,12,15,18,24,25,27,28,29,30,34,35,36,38,39,43,44,47,54,58,59,60,63,64,66,68,70,71,72,75,76,81,89,96,97,100,104,107,129,138,141,145,146,147,152,163,167,172,184,192,473,491,497,505
804,805,809,810,813,818,821,824,825,826,829,836,838,840,844,848,852,865,874

c = susie_rss(gwas$z, gwas$LD, n = unique(gwas$N), L = 10, max_iter = 2500, min_abs_corr = sqrt(0.4),verbose=TRUE,coverage = 0.99)
[1] "objective:-16969.693590593"
[1] "objective:-16969.6211774467"
[1] "objective:-16969.5931406855"
[1] "objective:-16969.5784032943"
[1] "objective:-16969.5687705576"
[1] "objective:-16969.5610090115"
[1] "objective:-16969.5536898941"
[1] "objective:-16969.5459284366"
[1] "objective:-16969.5364229281"
[1] "objective:-16969.5228911058"
[1] "objective:-16969.5017541661"
[1] "objective:-16969.4676286733"
[1] "objective:-16969.4152138044"
[1] "objective:-16969.3516179333"
[1] "objective:-16969.3016357743"
[1] "objective:-16969.2752213206"
[1] "objective:-16969.2639810721"
[1] "objective:-16969.2594363808"
[1] "objective:-16969.2573272553"
[1] "objective:-16969.2560825934"
[1] "objective:-16969.2552228556"
summary(c)

Variables in credible sets:

variable variable_prob cs
732 0.186933199 1
745 0.100456359 1
764 0.094084969 1
779 0.091830763 1
781 0.091798942 1
730 0.089566402 1
722 0.076941801 1
724 0.073764866 1
544 0.064415496 1
112 0.054626782 1
782 0.037255938 1
460 0.031753494 1
815 0.030612106 1
784 0.008439453 1
1350 0.007984235 1
1140 0.007260920 1
1197 0.006880256 1
1266 0.006870096 1
1315 0.006855586 1
712 0.006746137 1

Credible sets summary:

cs cs_log10bf cs_avg_r2 cs_min_r2
1 60.48482 0.9437577 0.8504178
variable
112,460,544,712,722,724,730,732,745,764,779,781,782,784,815,1140,1197,1266,1315,1350

chr1swallace commented
catheriz commented

str(gwas)
List of 12
$ snp : chr [1:1392] "chr6:31413077:G:T" "chr6:31413120:C:T" "chr6:31413166:G:A" "chr6:31413173:C:T" ...
$ chromosome: num [1:1392] 6 6 6 6 6 6 6 6 6 6 ...
$ position : num [1:1392] 31413077 31413120 31413166 31413173 31413197 ...
$ allele1 : chr [1:1392] "T" "T" "A" "T" ...
$ allele2 : chr [1:1392] "G" "C" "G" "C" ...
$ beta : num [1:1392] -0.354 -0.754 -0.358 -0.337 -0.354 ...
$ se : num [1:1392] 0.0942 0.1866 0.0941 0.1041 0.0941 ...
$ z : num [1:1392] -3.76 -4.04 -3.81 -3.23 -3.77 ...
$ N : int [1:1392] 12085 12085 12085 12085 12085 12085 12085 12085 12085 12085 ...
$ varbeta : num [1:1392] 0.00887 0.03481 0.00885 0.01084 0.00886 ...
$ LD : num [1:1392, 1:1392] 1 -0.154 0.999 0.841 1 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : chr [1:1392] "chr6:31413077:G:T" "chr6:31413120:C:T" "chr6:31413166:G:A" "chr6:31413173:C:T" ...
.. ..$ : chr [1:1392] "chr6:31413077:G:T" "chr6:31413120:C:T" "chr6:31413166:G:A" "chr6:31413173:C:T" ...
$ type : chr "cc"

laleoarrow commented

I've encountered a similar issue. It's possible that the susieR method did not achieve convergence, yet it could still produce a result. Conversely, when using the coloc-susie method, incorrect results may be yielded or can not be produced if the analysis is not properly converged？