add param min.maf to autoSVD

Author: privefl

DESCRIPTION

@@ -2,8 +2,8 @@ Encoding: UTF-8
 Package: bigsnpr
 Type: Package
 Title: Analysis of Massive SNP Arrays
-Version: 1.12.13
-Date: 2024-09-04
+Version: 1.12.14
+Date: 2024-09-10
 Authors@R: c(
     person("Florian", "Privé", email = "[email protected]", role = c("aut", "cre")),
     person("Michael", "Blum", role = "ths"),

NEWS.md

@@ -1,3 +1,7 @@
+## bigsnpr 1.12.14
+
+- Add a new `min.maf = 0.02` parameter to `snp_autoSVD()` and `bed_autoSVD()`. Then variants are now discarded when they have either a small MAC or a small MAF.
+
 ## bigsnpr 1.12.13
 
 - Can now use matrix accessors for class `bed_light` as well.

R/autoSVD.R

@@ -45,7 +45,9 @@ getIntervals <- function(x, n = 2) {
 #' @param alpha.tukey Default is `0.1`. The type-I error rate in outlier
 #'   detection (that is further corrected for multiple testing).
 #' @param min.mac Minimum minor allele count (MAC) for variants to be included.
-#'   Default is `10`.
+#'   Default is `10`. Can actually be higher because of `min.maf`.
+#' @param min.maf Minimum minor allele frequency (MAF) for variants to be included.
+#'   Default is `0.02`. Can actually be higher because of `min.mac`.
 #' @param max.iter Maximum number of iterations of outlier detection.
 #'   Default is `5`.
 #'
@@ -75,6 +77,7 @@ snp_autoSVD <- function(G,
                         int.min.size = 20,
                         alpha.tukey = 0.05,
                         min.mac = 10,
+                        min.maf = 0.02,
                         max.iter = 5,
                         is.size.in.bp = NULL,
                         ncores = 1,
@@ -90,15 +93,15 @@ snp_autoSVD <- function(G,
   # Verbose?
   printf2 <- function(...) if (verbose) printf(...)
 
-  if (min.mac > 0) {
+  if (min.mac > 0 && min.maf > 0) {
     maf <- snp_MAF(G, ind.row, ind.col, ncores = ncores)
-    min.maf <- min.mac / (2 * length(ind.row))
-    mac.nok <- (maf < min.maf)
-    printf2("Discarding %d variant%s with MAC < %s.\n", sum(mac.nok),
-            `if`(sum(mac.nok) > 1, "s", ""), min.mac)
-    ind.keep <- ind.col[!mac.nok]
+    maf.nok <- (maf < max(min.maf, min.mac / (2 * length(ind.row))))
+    printf2("Discarding %d variant%s with MAC < %s or MAF < %s.\n",
+            sum(maf.nok), `if`(sum(maf.nok) > 1, "s", ""), min.mac, min.maf)
+    ind.keep <- ind.col[!maf.nok]
   } else {
-    stop2("You cannot use variants with no variation; set min.mac > 0.")
+    stop2("You cannot use variants with no variation; %s",
+          "set min.mac > 0 and min.maf > 0.")
   }
 
   # First clumping
@@ -231,6 +234,7 @@ bed_autoSVD <- function(obj.bed,
                         int.min.size = 20,
                         alpha.tukey = 0.05,
                         min.mac = 10,
+                        min.maf = 0.02,
                         max.iter = 5,
                         ncores = 1,
                         verbose = TRUE) {
@@ -243,14 +247,15 @@ bed_autoSVD <- function(obj.bed,
   # Verbose?
   printf2 <- function(...) if (verbose) printf(...)
 
-  if (min.mac > 0) {
-    mac <- bed_MAF(obj.bed, ind.row, ind.col, ncores = ncores)$mac
-    mac.nok <- (mac < min.mac)
-    printf2("Discarding %d variant%s with MAC < %s.\n", sum(mac.nok),
-            `if`(sum(mac.nok) > 1, "s", ""), min.mac)
-    ind.keep <- ind.col[!mac.nok]
+  if (min.mac > 0 && min.maf > 0) {
+    info <- bed_MAF(obj.bed, ind.row, ind.col, ncores = ncores)
+    maf.nok <- (info$mac < min.mac | info$maf < min.maf)
+    printf2("Discarding %d variant%s with MAC < %s or MAF < %s.\n",
+            sum(maf.nok), `if`(sum(maf.nok) > 1, "s", ""), min.mac, min.maf)
+    ind.keep <- ind.col[!maf.nok]
   } else {
-    stop2("You cannot use variants with no variation; set min.mac > 0.")
+    stop2("You cannot use variants with no variation; %s",
+          "set min.mac > 0 and min.maf > 0.")
   }
 
   # First clumping

docs/news/index.html

@@ -10,6 +10,8 @@ CHR <- test$map$chromosome
 POS <- test$map$physical.pos / 10
 POS2 <- round(POS + 1)
 
+obj.bed <- bed(snp_writeBed(test, bedfile = tempfile(fileext = ".bed")))
+
 ################################################################################
 
 test_that("snp_autoSVD() works", {
@@ -65,8 +67,6 @@ test_that("snp_autoSVD() works", {
 
 test_that("bed_autoSVD() works", {
 
-  obj.bed <- bed(snp_writeBed(test, bedfile = tempfile(fileext = ".bed")))
-
   expect_error(bed_autoSVD(obj.bed, min.mac = 0), "no variation; set min.mac > 0")
   expect_output(bed_autoSVD(obj.bed, thr.r2 = NA), "Skipping clumping.")
 
@@ -96,3 +96,28 @@ test_that("bed_autoSVD() works", {
 })
 
 ################################################################################
+
+test_that("MAC/MAF thresholds work", {
+
+  info <- bed_MAF(obj.bed)
+
+  replicate(10, {
+
+    min.mac <- sample(1:40, size = 1)
+    min.maf <- runif(n = 1, min = 0.01, max = 0.1)
+
+    obj.svd1 <- snp_autoSVD(G, CHR, size = 5, min.mac = min.mac, min.maf = min.maf,
+                            thr.r2 = NA, max.iter = 0, verbose = FALSE)
+    ind1 <- attr(obj.svd1, "subset")
+    expect_true(all(info$maf[ind1] >= min.maf))
+    expect_true(all(info$mac[ind1] >= min.mac))
+
+    obj.svd2 <- bed_autoSVD(obj.bed, min.mac = min.mac, min.maf = min.maf,
+                            thr.r2 = NA, max.iter = 0, verbose = FALSE)
+    ind2 <- attr(obj.svd2, "subset")
+    expect_true(all(info$maf[ind2] >= min.maf))
+    expect_true(all(info$mac[ind2] >= min.mac))
+  })
+})
+
+################################################################################