coloc icon indicating copy to clipboard operation
coloc copied to clipboard

Published 20 hours ago • verified publisher icon chr1swallace

runsusie and susieR results differ

Open maguileraf opened this issue 1 year ago • 8 comments

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?

maguileraf avatar Nov 28 '23 00:11 maguileraf

No, this shouldn't happen because runsusie calls susieR. Could you please share the data for this?

-- https://chr1swallace.github.io

From: Marcela Johnson @.> Sent: Tuesday, November 28, 2023 12:55:42 AM To: chr1swallace/coloc @.> Cc: Subscribed @.***> Subject: [chr1swallace/coloc] runsusie and susieR results differ (Issue #142)

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?

— Reply to this email directly, view it on GitHubhttps://github.com/chr1swallace/coloc/issues/142, or unsubscribehttps://github.com/notifications/unsubscribe-auth/AAQWR2ERRKI2WCX3NJD3N23YGUZA5AVCNFSM6AAAAAA744DHHOVHI2DSMVQWIX3LMV43ASLTON2WKOZSGAYTGMZZGM4DANY. You are receiving this because you are subscribed to this thread.Message ID: @.***>

chr1swallace avatar Nov 28 '23 07:11 chr1swallace

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

maguileraf avatar Nov 28 '23 11:11 maguileraf

The z scores in runsusie are beta/sqrt(varbeta)

What is the magnitude of your z scores at those snps though? It is very unusual to fine 6 single variant credible sets

From: Marcela Johnson @.> Sent: Tuesday, November 28, 2023 11:25:45 AM To: chr1swallace/coloc @.> Cc: Chris Wallace @.>; Comment @.> Subject: Re: [chr1swallace/coloc] runsusie and susieR results differ (Issue #142)

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) They are both using the same LD matrix and same data. However, maybe the Z scores are being calculated differently?

Screenshot.2023-11-28.at.6.24.08.AM.png (view on web)https://github.com/chr1swallace/coloc/assets/49694113/389dc460-ede4-424d-8ad9-173c0fc46e21 Screenshot.2023-11-28.at.6.24.26.AM.png (view on web)https://github.com/chr1swallace/coloc/assets/49694113/e886dd55-8312-466f-bfe3-034b3603a215 , refine=T)

— Reply to this email directly, view it on GitHubhttps://github.com/chr1swallace/coloc/issues/142#issuecomment-1829631369, or unsubscribehttps://github.com/notifications/unsubscribe-auth/AAQWR2ANRSNNKQLYDEDIZG3YGXC3TAVCNFSM6AAAAAA744DHHOVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMYTQMRZGYZTCMZWHE. You are receiving this because you commented.

chr1swallace avatar Nov 28 '23 11:11 chr1swallace

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

maguileraf avatar Nov 28 '23 13:11 maguileraf

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

catheriz avatar Mar 07 '24 20:03 catheriz

Can you show me str(gwas)?

-- https://chr1swallace.github.io

From: catheriz @.> Sent: Thursday, March 7, 2024 8:39:46 PM To: chr1swallace/coloc @.> Cc: Chris Wallace @.>; Comment @.> Subject: Re: [chr1swallace/coloc] runsusie and susieR results differ (Issue #142)

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

— Reply to this email directly, view it on GitHubhttps://github.com/chr1swallace/coloc/issues/142#issuecomment-1984389925, or unsubscribehttps://github.com/notifications/unsubscribe-auth/AAQWR2GUF5FI4VMVKYMJUMLYXDGBFAVCNFSM6AAAAAA744DHHOVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMYTSOBUGM4DSOJSGU. You are receiving this because you commented.Message ID: @.***>

chr1swallace avatar Mar 07 '24 21:03 chr1swallace

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"

catheriz avatar Mar 07 '24 21:03 catheriz

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？

laleoarrow avatar Apr 15 '24 01:04 laleoarrow