myles-lewis / locuszoomr

A pure R implementation of locuszoom for plotting genetic data at genomic loci accompanied by gene annotations.
GNU General Public License v3.0
18 stars 5 forks source link

link_LD matches 0 SNPs #23

Closed msbrasher closed 2 months ago

msbrasher commented 2 months ago

Hi Myles

I am trying to create plots for loci on the X-chromosome, and I'm running into trouble with adding LD information where the link_LD() function is running, but returning zero matched SNPs (see output below).

Screenshot 2024-07-12 at 1 49 33 PM

The output plots look great other than the missing LD information: Screenshot 2024-07-12 at 1 51 56 PM

Here is the code to reproduce (and test file gwas.txt):

library(locuszoomr)
library(data.table)
library(EnsDb.Hsapiens.v86)

test_gwas <- fread("gwas.txt")

token <- my_token

pks <- quick_peak(test_gwas, chrom = "CHROM", pos = "GENPOS", p = "P")
top_snps <- test_gwas$rsid[pks]
head(top_snps)

all_loci <- lapply(top_snps, function(i) {
  loc <- locus(data = test_gwas, index_snp = i, fix_window = 1e6, p = "P",
               ens_db = "EnsDb.Hsapiens.v86")

  loc <- link_LD(loc, pop = "EUR", token = token, genome_build = "grch38_high_coverage")

  link_recomb(loc, genome = "hg38")
})

pdf("sle_loci.pdf")
tmp <- lapply(all_loci, locus_plot, labels = "index")
dev.off()

The comparable query directly to the LDproxy function in the LDlinkR package does return LD results out <- LDproxy(snp = "chrX:55627300", pop = "EUR", token = token, genome_build = "grch38_high_coverage")

For reference: The GWAS summary stats I have are from TOPMed imputed data, so they are hg38, and the index SNPs I'm interested in are only present in the grch38_high_coverage build. I'm not sure if using that build might be causing the problem? I'm using locuszoomr version 0.3.1

myles-lewis commented 2 months ago

Hi @msbrasher,

Thanks for raising this. I've pushed a fix to the code for link_LD(). It should be able to cope with your example now. Note this only seems to work for the LDproxy method. LDlinkR::LDmatrix() returns rs IDs even if you provide chromosome coordinate SNP ids. So matching will fail within the locus object. Let me know if it works ok.

Bw, Myles

msbrasher commented 2 months ago

Hi Myles,

Looks like it's working great now! Thank you so much for the fast response!

Maizy