fix(gwas): make adjustPval and manhatan plot handle missing p values

Author: juliendiot42

src/gwas.R

@@ -146,7 +146,7 @@ adjustPval <- function(p, adj_method, thresh_p = NULL) {
 
   # check p
   logger$log("Check p values ...")
-  if (any(p < 0) || any(p > 1)) {
+  if (any(p < 0, na.rm = TRUE) || any(p > 1, na.rm = TRUE)) {
     # The `p.adjust` function will proceed to the calculation even if the p-values
     # are invalid. To avoid unexpected behaviour, I'd prefer this function
     # to crash in such case.

src/plot.R

@@ -37,6 +37,7 @@ manPlot <- function(gwas,
     all(is.na(x))
   })
   gwas <- gwas[!allNaLines, ]
+  gwas <- gwas[!is.na(gwas$p), ]
   logger$log("Remove NAs DONE")
 
   # P-Values adjustment ----
@@ -81,7 +82,7 @@ manPlot <- function(gwas,
   # get significant SNP ----
   if (highlightSinif) {
     logger$log("Extract significant SNP ...")
-    significantSNP <- gwas[gwas$p_adj <= thresh_p, "id"]
+    significantSNP <- gwas[ifelse(is.na(gwas$p_adj), FALSE, gwas$p_adj <= thresh_p), "id"]
     if (length(significantSNP) == 0) {
       significantSNP <- NULL
     }

tests/testthat/test_2_gwas.R

@@ -155,11 +155,36 @@ capture.output({
   }
 
   # adjustPval with edge case
-  test_that("adjustPval handle gwas with missing p values", {
+  test_that("adjustPval handle gwas with some missing p values", {
     n_p <- 100
     p <- runif(n_p)
-    p[sample.int(n_p, 5)] <- runif(5, min = 0, max = 0.05)
-    p[sample.int(n_p, 5)] <- NA
+
+    missing_pval_indices <- sample.int(n_p, 5)
+    p[missing_pval_indices] <- NA
+    thresh_p <- 0.05
+
+    for (adjMeth in c(
+      "holm",
+      "hochberg",
+      "bonferroni",
+      "BH",
+      "BY",
+      "fdr",
+      "none"
+    )) {
+      expect_no_error({
+        adj <- adjustPval(
+          p = p,
+          adj_method = adjMeth,
+          thresh_p = thresh_p
+        )
+      })
+      expect_true(all(is.na(adj$p_adj[missing_pval_indices])))
+    }
+  })
+
+  test_that("adjustPval handle gwas with all missing p values", {
+    p <- rep(NA, 100)
     thresh_p <- 0.05
 
     for (adjMeth in c(
@@ -178,6 +203,7 @@ capture.output({
           thresh_p = thresh_p
         )
       })
+      expect_true(all(is.na(adj$p_adj)))
     }
   })
 })

tests/testthat/test_4_plots.R

@@ -71,6 +71,53 @@ capture_output({
                 expect_true(is.null(mP))
               }
             })
+
+            test_that(paste(testName, "with some missing p values"), {
+              resGwas$p[c(2, 4, 10)] <- NA
+              expect_no_error({
+                mP <- manPlot(
+                  gwas = resGwas,
+                  adj_method = "bonferroni",
+                  thresh_p = 0.05,
+                  chr = chr,
+                  interactive = interactive,
+                  filter_pAdj = 1,
+                  filter_nPoints = Inf,
+                  filter_quant = 1,
+                  title = "test manPlot"
+                )
+              })
+              if (interactive) {
+                expect_true(all.equal(class(mP), c("plotly", "htmlwidget")))
+              } else {
+                expect_true(is.null(mP))
+              }
+            })
+
+            test_that(paste(testName, "with all missing p values"), {
+              resGwas$p <- NA
+              expect_warning(
+                {
+                  mP <- manPlot(
+                    gwas = resGwas,
+                    adj_method = "bonferroni",
+                    thresh_p = 0.05,
+                    chr = chr,
+                    interactive = interactive,
+                    filter_pAdj = 1,
+                    filter_nPoints = Inf,
+                    filter_quant = 1,
+                    title = "test manPlot"
+                  )
+                },
+                regexp = "There is no points to display"
+              )
+              if (interactive) {
+                expect_true(all.equal(class(mP), c("plotly", "htmlwidget")))
+              } else {
+                expect_true(is.null(mP))
+              }
+            })
           }
         }
       }