search

privefl / paper-ldpred2

Paper discribing LDpred2
https://doi.org/10.1093/bioinformatics/btaa1029
16 stars 12 forks source link

Issue with example-with-provided-ldref.R #14

Closed lm380-r closed 1 year ago

lm380-r commented 1 year ago

Hi,

First of all, thank you for this tool. I was trying to work with the example-with-provided-ldref.R script but when I keep getting an error when I try to run lines 50-70.

Code: tmp <- tempfile(tmpdir = "tmp-data")

for (chr in 1:22) {

cat(chr, ".. ", sep = "")

ind.chr <- which(df_beta$chr == chr) ind.chr2 <- df_beta$_NUM_ID_[ind.chr] ind.chr3 <- match(ind.chr2, which(map_ldref$chr == chr))

corr_chr <- readRDS(paste0("ld_chr", chr, ".rds"))[ind.chr3, ind.chr3]

if (chr == 1) { corr <- as_SFBM(corr_chr, tmp, compact = TRUE) } else { corr$add_columns(corr_chr, nrow(corr)) } }

Error: 1.. Error in gzfile(file, "rb") : cannot open the connection In addition: Warning message: In gzfile(file, "rb") : cannot open compressed file 'ld_chr1.rds', probable reason 'No such file or directory'

I am new to this so any help would be greatly appreciated. I can get ind.chr, ind.chr2 and ind.chr3 commands to run but I'm not sure how to save those so I can use them for the rest of the code.

Thank you in advance!

privefl commented 1 year ago

You need to replace with the path where you downloaded the LD matrices.

lm380-r commented 1 year ago

Thank you!

One last question (sorry!).. any idea what this means?

1.. Error in intI(i, n = d[1L], dn[[1L]], give.dn = FALSE) : index larger than maximal 87145

I tried to google it... I know the error comes from this line:

corr_chr <- readRDS(paste0("LD_chr", chr, ".rds"))[ind.chr3, ind.chr3]

privefl commented 1 year ago

You have indices in ind.chr3 that are larger to the dimensions of the matrix you get from readRDS(paste0("LD_chr", chr, ".rds")).

lm380-r commented 1 year ago

The issue was because I was using map_hm3_plus.rds instead of map.rds. The most recent tutorial recommended using map_hm3_plus.rds so I thought I should use that but I learned I cannot.

Thank you!

lm380-r commented 1 year ago

Sorry to follow up, but where would I find the LD block files (by chromosome) for this step:

corr_chr <- readRDS(paste0("data/corr_hm3_plus/LD_with_blocks_chr", chr, ".rds"))[ind.chr3, ind.chr3]

For the hapmap+ SNPs?

The tutorial only points to the map_hm3_plus.rds but the code requires LD blocks by chromosome and if you try to use the precomputed LD blocks based on the previous hapmap SNPs, the code errors out because there are more SNPs provided in the .rds file than in the LD blocks.

privefl commented 1 year ago

In figshare, there is a G-drive link in the description.