coloc.abf returns NaN for all posteriors

[Bonder-MJ]

Hello,

Thanks for the great software! I was trying to use it on my data and it seems that I encounterd a situation where all posteriors are returned NaN but there is no reason / error reported for this.

I checked the "check_dataset()" outputs for both the inputs, both report "## NULL" and I do get SNP priors & hypothesis priors but no posteriors.

What I'm trying to run:

`> coloc.abf(dataset1=qtlD, dataset2=GwasCC)
PP.H0.abf PP.H1.abf PP.H2.abf PP.H3.abf PP.H4.abf
NaN NaN NaN NaN NaN
[1] "PP abf for shared variant: NaN%"
Coloc analysis of trait 1, trait 2

SNP Priors
p1 p2 p12
1e-04 1e-04 1e-05

Hypothesis Priors
H0 H1 H2 H3 H4
0.7534289 0.1115 0.1115 0.0124211 0.01115

Posterior
nsnps H0 H1 H2 H3 H4
1115 NaN NaN NaN NaN NaN `

I'm not sure what is going on here and I hoped you might have a suggestion!

Thanks,
Marc

[chr1swallace]

can you show me str() on each dataset?

…

On Mon, 2021-11-29 at 07:19 -0800, Marc Jan Bonder wrote: Hello, Thanks for the great software! I was trying to use it on my data and it seems that I encounterd a situation where all posteriors are returned NaN but there is no reason / error reported for this. I checked the "check_dataset()" outputs for both the inputs, both report "## NULL" and I do get SNP priors & hypothesis priors but no posteriors. What I'm trying to run: `> coloc.abf(dataset1=qtlD, dataset2=GwasCC) PP.H0.abf PP.H1.abf PP.H2.abf PP.H3.abf PP.H4.abf NaN NaN NaN NaN NaN [1] "PP abf for shared variant: NaN%" Coloc analysis of trait 1, trait 2 SNP Priors p1 p2 p12 1e-04 1e-04 1e-05 Hypothesis Priors H0 H1 H2 H3 H4 0.7534289 0.1115 0.1115 0.0124211 0.01115 Posterior nsnps H0 H1 H2 H3 H4 1115 NaN NaN NaN NaN NaN ` I'm not sure what is going on here and I hoped you might have a suggestion! Thanks, Marc — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub, or unsubscribe. Triage notifications on the go with GitHub Mobile for iOS or Android.

[Bonder-MJ]

Thanks for the super fast reply. Here you go:

str(GWASStats)
List of 8
$ beta : Named num [1:1115] 1.225 -0.434 1.748 1.748 -1.386 ...
..- attr(, "names")= chr [1:1115] "6:32385399_G_A" "6:32385453_G_A" "6:32385470_G_A" "6:32385573_A_G" ...
$ varbeta : Named num [1:1115] 0.0829 0.1545 0.2113 0.2113 0.2874 ...
..- attr(, "names")= chr [1:1115] "6:32385399_G_A" "6:32385453_G_A" "6:32385470_G_A" "6:32385573_A_G" ...
$ snp : chr [1:1115] "6:32385399_G_A" "6:32385453_G_A" "6:32385470_G_A" "6:32385573_A_G" ...
$ position: int [1:1115] 32385399 32385453 32385470 32385573 32385596 32385605 32385653 32385750 32385782 32385873 ...
$ type : chr "quant"
$ sdY : num 1
$ N : int [1:1115] 42 42 42 42 42 42 42 42 42 42 ...
$ MAF : num [1:1115] 0.298 0.286 0.167 0.167 0.119 ...

str(geneStats)
List of 5
$ beta : Named num [1:1115] -0.779 0.447 -0.144 -0.144 0.32 ...
..- attr(, "names")= chr [1:1115] "6:32385399_G_A" "6:32385453_G_A" "6:32385470_G_A" "6:32385573_A_G" ...
$ varbeta : Named num [1:1115] 0 0.000324 0.000504 0.000504 0.000833 ...
..- attr(, "names")= chr [1:1115] "6:32385399_G_A" "6:32385453_G_A" "6:32385470_G_A" "6:32385573_A_G" ...
$ snp : chr [1:1115] "6:32385399_G_A" "6:32385453_G_A" "6:32385470_G_A" "6:32385573_A_G" ...
$ position: int [1:1115] 32385399 32385453 32385470 32385573 32385596 32385605 32385653 32385750 32385782 32385873 ...
$ type : chr "cc"`

And here is the combined 'plot_dataset" plot

As you can see it is a locus on Chr 6 so it might have to do with complex LD around the HLA?

[chr1swallace]

perhaps it is the extreme p values, looking in the HLA region? to check, could you do finemap.abf() on each dataset? If only the one with tiny p values gives NA, then that would support that hypothesis, but I am struggling looking at the code to see how that would cause errors. Two things. 1. I want to fix this, coloc should not return nan without warnings first. Could you possibly share your data? 2. I would be skeptical about any analysis that invokes the single causal variant assumption in the HLA region.

…

On Mon, 2021-11-29 at 07:31 -0800, Marc Jan Bonder wrote: Thanks for the super fast reply. Here you go: `> str(GWASStats) 'data.frame': 1115 obs. of 13 variables: $ CHR : int 6 6 6 6 6 6 6 6 6 6 ... $ BP : int 32385399 32385453 32385470 32385573 32385596 32385605 32385653 32385750 32385782 32385873 ... $ SNP : chr "rs3135377" "rs9268544" "rs3135376" "rs3135375" ... $ A1 : chr "G" "G" "G" "A" ... $ A2 : chr "A" "A" "A" "G" ... $ N : int 15 15 15 15 15 15 15 15 15 15 ... $ OR : num 0.459 1.563 0.866 0.866 1.378 ... $ P : num 0.00 4.21e-136 1.50e-10 1.54e-10 1.16e-28 ... $ OR.R. : num 0.459 1.56 0.872 0.872 1.379 ... $ P.R. : num 2.30e-186 1.23e-39 1.33e-03 1.33e-03 4.80e-27 ... $ Q : num 0.0353 0.0001 0 0 0.4013 ... $ I : num 43.83 68.36 69.36 69.31 4.55 ... $ markerId: chr "6:32385399_G_A" "6:32385453_G_A" "6:32385470_G_A" "6:32385573_A_G" ... > str(geneStats) > 'data.frame': 1115 obs. of 29 variables: > $ feature_id : chr "HLA-DRB5" "HLA-DRB5" "HLA-DRB5" "HLA-DRB5" ... > $ snp_id : chr "6:32385399" "6:32385453" "6:32385470" "6:32385573" > ... > $ p_value : num 2.09e-05 2.70e-01 1.44e-04 1.44e-04 9.71e-03 ... > $ beta : num 1.225 -0.434 1.748 1.748 -1.386 ... > $ beta_se : num 0.288 0.393 0.46 0.46 0.536 ... > $ empirical_feature_p_value: num 0.0134 1 0.0435 0.0435 0.5078 ... > $ feature_chromosome : int 6 6 6 6 6 6 6 6 6 6 ... > $ feature_start : int 32485120 32485120 32485120 32485120 32485120 > 32485120 32485120 32485120 32485120 32485120 ... > $ feature_end : int 32498064 32498064 32498064 32498064 32498064 > 32498064 32498064 32498064 32498064 32498064 ... > $ strand : int -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 ... > $ Gene.name : chr "HLA-DRB5" "HLA-DRB5" "HLA-DRB5" "HLA-DRB5" ... > $ Gene.stable.ID : chr "ENSG00000198502" "ENSG00000198502" > "ENSG00000198502" "ENSG00000198502" ... > $ Gene.type : chr "protein_coding" "protein_coding" > "protein_coding" "protein_coding" ... > $ duplicted : int 0 0 0 0 0 0 0 0 0 0 ... > $ n_samples : int 42 42 42 42 42 42 42 42 42 42 ... > $ n_e_samples : int 42 42 42 42 42 42 42 42 42 42 ... > $ alpha_param : num 0.611 0.611 0.611 0.611 0.611 ... > $ beta_param : num 34.6 34.6 34.6 34.6 34.6 ... > $ rho : num 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 ... > $ snp_chromosome : int 6 6 6 6 6 6 6 6 6 6 ... > $ snp_position : int 32385399 32385453 32385470 32385573 32385596 > 32385605 32385653 32385750 32385782 32385873 ... > $ assessed_allele : chr "G" "A" "G" "A" ... > $ call_rate : num 1 1 1 1 1 1 1 1 1 1 ... > $ maf : num 0.298 0.286 0.167 0.167 0.119 ... > $ hwe_p : num 0.133 0.454 0.572 0.572 1 ... > $ QTL : chr "6:32385399-HLA-DRB5" "6:32385453-HLA-DRB5" > "6:32385470-HLA-DRB5" "6:32385573-HLA-DRB5" ... > $ A1 : chr "G" "G" "G" "A" ... > $ A2 : chr "A" "A" "A" "G" ... > $ markerId : chr "6:32385399_G_A" "6:32385453_G_A" "6:32385470_G_A" > "6:32385573_A_G" ...` And here is the combined 'plot_dataset" plot — You are receiving this because you commented. Reply to this email directly, view it on GitHub, or unsubscribe. Triage notifications on the go with GitHub Mobile for iOS or Android.

[Bonder-MJ]

Hi,

Sorry for only getting back to you now. I when back and double checked all input data again and I found that there are entries with P of 0 and beta_se of 0, after filtering out these unexpected entries in my datasets it works all fine!

Thanks for your quick replies and help!

Best,
Marc

[chr1swallace]

thanks for getting back to me, I can add a check for 0 (or negative) betavar.

…

On Wed, 2021-12-01 at 12:48 -0800, Marc Jan Bonder wrote: Hi, Sorry for only getting back to you now. I when back and double checked all input data again and I found that there are entries with P of 0 and beta_se of 0, after filtering out these unexpected entries in my datasets it works all fine! Thanks for your quick replies and help! Best, Marc — You are receiving this because you commented. Reply to this email directly, view it on GitHub, or unsubscribe. Triage notifications on the go with GitHub Mobile for iOS or Android.