Add parameter rsid to snp_asGeneticPos()

DESCRIPTION

@@ -2,8 +2,8 @@ Encoding: UTF-8
 Package: bigsnpr
 Type: Package
 Title: Analysis of Massive SNP Arrays
-Version: 1.7.0
-Date: 2021-03-16
+Version: 1.7.1
+Date: 2021-03-30
 Authors@R: c(
     person("Florian", "Privé", 
            email = "[email protected]", 

NEWS.md

@@ -1,3 +1,11 @@
+## bigsnpr 1.7.1
+
+- Add parameter `rsid` to `snp_asGeneticPos()` to also allow matching with rsIDs.
+
+## bigsnpr 1.7.0
+
+- Add function `snp_lassosum2()` to train the lassosum models using the exact same input data as LDpred2.
+
 ## bigsnpr 1.6.7
 
 - Add parameter `report_step` in `snp_ldpred2_auto()` to report some of the internal sampling betas.

R/match-alleles.R

@@ -260,16 +260,22 @@ same_ref <- function(ref1, alt1, ref2, alt2) {
 #' @inheritParams bigsnpr-package
 #' @param dir Directory where to download and decompress files.
 #'   Default is `tempdir()`. Directly use files there if already present.
+#' @param rsid If providing rsIDs, the matching is performed using those
+#'   (instead of positions) and variants not matched are interpolated using
+#'   spline interpolation of variants that have been matched.
 #'
 #' @return The new vector of genetic positions.
 #' @export
 #'
-snp_asGeneticPos <- function(infos.chr, infos.pos, dir = tempdir(), ncores = 1) {
+snp_asGeneticPos <- function(infos.chr, infos.pos, dir = tempdir(), ncores = 1,
+                             rsid = NULL) {
 
   assert_package("R.utils")
   assert_lengths(infos.chr, infos.pos)
+  if (!is.null(rsid)) assert_lengths(rsid, infos.pos)
+
+  snp_split(infos.chr, function(ind.chr, pos, dir, rsid) {
 
-  snp_split(infos.chr, function(ind.chr, pos, dir) {
     chr <- attr(ind.chr, "chr")
     basename <- paste0("chr", chr, ".OMNI.interpolated_genetic_map")
     mapfile <- file.path(dir, basename)
@@ -281,10 +287,26 @@ snp_asGeneticPos <- function(infos.chr, infos.pos, dir = tempdir(), ncores = 1)
       R.utils::gunzip(gzfile)
     }
     map.chr <- bigreadr::fread2(mapfile, showProgress = FALSE)
-    ind <- bigutilsr::knn_parallel(as.matrix(map.chr$V2), as.matrix(pos[ind.chr]),
-                                   k = 1, ncores = 1)$nn.idx
-    map.chr$V3[ind]
-  }, combine = "c", pos = infos.pos, dir = dir, ncores = ncores)
+
+    if (is.null(rsid)) {
+      ind <- bigutilsr::knn_parallel(as.matrix(map.chr$V2), as.matrix(pos[ind.chr]),
+                                     k = 1, ncores = 1)$nn.idx
+      new_pos <- map.chr$V3[ind]
+    } else {
+      ind <- match(rsid[ind.chr], map.chr$V1)
+      new_pos <- map.chr$V3[ind]
+
+      indNA <- which(is.na(ind))
+      if (length(indNA) > 0) {
+        pos.chr <- pos[ind.chr]
+        new_pos[indNA] <- suppressWarnings(
+          stats::spline(pos.chr, new_pos, xout = pos.chr[indNA], method = "hyman")$y)
+      }
+    }
+
+    new_pos
+
+  }, combine = "c", pos = infos.pos, dir = dir, rsid = rsid, ncores = ncores)
 }
 
 ################################################################################

tests/testthat/test-5-match.R

@@ -69,3 +69,29 @@ expect_equal(same_ref(ref1 = info_snp$a1, alt1 = info_snp$a0,
              c(TRUE, FALSE, TRUE, TRUE, TRUE))
 
 ################################################################################
+
+
+test_that("snp_asGeneticPos() works", {
+
+  info <- data.frame(
+    chr = rep(1L, 5L),
+    pos = c(853954L, 854250L, 864938L, 870645L, 873558L),
+    rsid = c("rs1806509", "rs7537756", "rs2340587", "rs28576697", "rs1110052"))
+
+  map <- data.frame(
+    V1 = c("rs1806509", "rs7537756", "rs28576697", "rs1110052"),
+    V2 = c(853954L, 854250L, 870645L, 873558L),
+    V3 = c(0.194323402834, 0.194576977815, 0.202835640491, 0.203874368612))
+  bigreadr::fwrite2(map, file.path(tempdir(), "chr1.OMNI.interpolated_genetic_map"))
+
+  res1 <- snp_asGeneticPos(info$chr, info$pos, dir = tempdir(), ncores = 2)
+  expect_equal(res1, map$V3[c(1, 2, 3, 3, 4)])
+
+  res2 <- snp_asGeneticPos(info$chr, info$pos, dir = tempdir(), ncores = 2,
+                           rsid = info$rsid)
+  expect_equal(res2[-3], map$V3)
+  expect_gt(res2[3], res2[2])
+  expect_lt(res2[3], res2[4])
+})
+
+################################################################################