chr1swallace / coloc

Repo for the R package coloc
144 stars 44 forks source link

Update checks #62

Closed jdblischak closed 2 years ago

jdblischak commented 2 years ago

I was attempting to do some quick tests with runsusie(), but I couldn't get it to run despite it passing check_dataset() and check_ld(). The code below attempts to run the example from one of the susieR vignettes through runsusie() (so I could compare the results between runsusie() and susie_rss()).

library(coloc)
library(susieR)

data("N3finemapping")
attach(N3finemapping)
sumstats <- univariate_regression(X, Y[, 1])
z_scores <- sumstats$betahat / sumstats$sebetahat
R <- cor(X)
fitted_rss <- susie_rss(z_scores, R, L = 10)

N3coloc <- list(
  beta = sumstats$betahat,
  varbeta = sumstats$sebetahat^2,
  N = 574,
  type = "quant",
  LD = R,
  sdY = sd(Y[, 1])
)
check_dataset(N3coloc)
## NULL
coloc:::check_ld(N3coloc, R)
N3colocsusie <- runsusie(N3coloc)
## running max iterations: 100
## Error in (function (z, R, z_ld_weight = 0, L = 10, prior_variance = 50,  :
##   The dimension of correlation matrix (0 by 0) does not agree with expected (1001 by 1001)

My input dataset didn't have the required element "snp". This wasn't caught by check_dataset() because it only loops through nd, which is the names of the input list.

https://github.com/chr1swallace/coloc/blob/f0c4a6a8126cf95c254d81246e2fa74912131111/R/check.R#L71-L72

My LD matrix didn't have row or column names, but this wasn't caught because NULL is identical to NULL:

https://github.com/chr1swallace/coloc/blob/f0c4a6a8126cf95c254d81246e2fa74912131111/R/check.R#L162-L165

And the error from susie_rss() is because runsusie() filters/order the LD matrix to match the "snp" vector:

https://github.com/chr1swallace/coloc/blob/f0c4a6a8126cf95c254d81246e2fa74912131111/R/susie.R#L460

Which I suspect is the problem in Issue #59 (they have an element "SNP" but not "snp").

This PR updates check_dataset() and check_ld() to catch these errors before running susie_rss(). Please let me know if you'd like me to make any changes. Thanks for developing such a useful package!