handle log10 p

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: MungeSumstats
 Type: Package
 Title: Standardise summary statistics from GWAS
-Version: 1.13.1
+Version: 1.13.2
 Authors@R:
     c(person(given = "Alan",
            family = "Murphy",

NEWS.md

@@ -1,3 +1,8 @@
+## CHANGES IN VERSION 1.13.2
+
+### New features
+* Handling of -log10 p-values (outside of VCFs) added.
+
 ## CHANGES IN VERSION 1.13.1
 
 ### New features

R/format_sumstats.R

@@ -482,6 +482,10 @@ format_sumstats <- function(path,
                                 "Corrected"=rep("ES",14))
           mapping_file <- rbind(mapping_file,es_cols)
         }
+
+        #### Check 40:Check for log10 p instead of p ####
+        sumstats_return <-
+          read_log_pval(sumstats_dt = sumstats_return$sumstats_dt)
 
         #### Check 2:Check for effect direction ####
         sumstats_return <-

R/read_log_pval.R

@@ -0,0 +1,72 @@
+#' Read -log10 p-value column
+#' 
+#' Parse p-value column in VCF file.of other general -loq10 p-values
+#' 
+#' @param sumstats_dt Summary stats data.table. 
+#' @param return_list Binary, whether to return the dt in a list or not - list 
+#' is standard for the `format_sumstats()` function.
+#' @inheritParams format_sumstats 
+#' 
+#' @returns Null output.
+#' @keywords internal
+#' @importFrom data.table setnames
+read_log_pval <- function(sumstats_dt,mapping_file = sumstatsColHeaders,
+                          return_list=TRUE){
+    P <- NULL
+    # Need to convert P-value, currently -log10
+    #cgreate a check based on mapping file
+    column_headers <- names(sumstats_dt) 
+    column_headers_upper <- toupper(names(sumstats_dt))
+    #get all p mappings 
+    mapping_file_p <- mapping_file[mapping_file$Corrected=="P",]
+    #now create all possible -log10 p-value
+    log_p_col <- c(paste0("-LOG10_",mapping_file_p$Uncorrected),
+                   paste0("LOG10_",mapping_file_p$Uncorrected),
+                   paste0("-LOG10 ",mapping_file_p$Uncorrected),
+                   paste0("LOG10 ",mapping_file_p$Uncorrected),
+                   paste0("-LOG10-",mapping_file_p$Uncorrected),
+                   paste0("LOG10-",mapping_file_p$Uncorrected),
+                   paste0("-LOG10.",mapping_file_p$Uncorrected),
+                   paste0("LOG10.",mapping_file_p$Uncorrected),
+                   paste0("-LOG10",mapping_file_p$Uncorrected),
+                   paste0("LOG10",mapping_file_p$Uncorrected),
+                   paste0("LOG_",mapping_file_p$Uncorrected),
+                   paste0("-L_",mapping_file_p$Uncorrected),
+                   paste0("L_",mapping_file_p$Uncorrected),
+                   paste0("LOG-",mapping_file_p$Uncorrected),
+                   paste0("-L-",mapping_file_p$Uncorrected),
+                   paste0("L-",mapping_file_p$Uncorrected),
+                   paste0("LOG.",mapping_file_p$Uncorrected),
+                   paste0("-L.",mapping_file_p$Uncorrected),
+                   paste0("L.",mapping_file_p$Uncorrected),
+                   paste0("-L",mapping_file_p$Uncorrected),
+                   paste0("L",mapping_file_p$Uncorrected),
+                   "LP")
+    lp_found <- FALSE
+    for(col_i in  column_headers_upper){
+      if(isTRUE(any(col_i==log_p_col))){
+        lp_found<-TRUE
+        lp_col_up<-col_i
+      }
+    }
+    if (isTRUE(lp_found)) {
+        messager("sumstats has -log10 P-values; these will be",
+                 "converted to unadjusted p-values in the 'P' column.")
+        lp_col <- column_headers[[which(column_headers_upper %in% lp_col_up)]]
+        #check if -log10 or log10
+        if(max(sumstats_dt[,get(lp_col)])<=0){
+          #log10 - less likely so more lenient check
+          sumstats_dt[, P := 10^(as.numeric(get(lp_col)))]
+        }else{
+          #-log10 - more likely
+          sumstats_dt[, P := 10^(-1 * as.numeric(get(lp_col)))]  
+        }
+
+    }
+    #### Return format ####
+    if(return_list){
+      return(list("sumstats_dt" = sumstats_dt))
+    }else {
+      return(sumstats_dt)
+    }
+}

R/read_vcf.R

@@ -136,7 +136,8 @@ read_vcf <- function(path,
     #### Prepare SNP column ####
     read_vcf_markername(sumstats_dt = vcf_dt)
     #### Prepare/un-log P col ####
-    read_vcf_pval(sumstats_dt = vcf_dt)
+    sumstats_dt <- read_log_pval(sumstats_dt = vcf_dt,
+                                 return_list=FALSE)
     #### Prepare INFO col ####
     read_vcf_info(sumstats_dt = vcf_dt)
     #### Rename start col #####

tests/testthat/test-log10_p.R

@@ -0,0 +1,51 @@
+test_that("log10 p-value", {
+    ## Call uses reference genome as default with more than 2GB of memory,
+    ## which is more than what 32-bit Windows can handle so remove tests
+    is_32bit_windows <-
+        .Platform$OS.type == "windows" && .Platform$r_arch == "i386"
+    if (!is_32bit_windows) {
+      # read it in and make N
+      sumstats_dt <- data.table::fread(system.file("extdata", "eduAttainOkbay.txt",
+                                                   package = "MungeSumstats"), 
+                                       nThread = 1)
+      #test -log10
+      sumstats_dt[,LP:=-1*log(Pval,base=10)]
+      data.table::setnames(sumstats_dt,"Pval","Pval_org")
+      reformatted <- MungeSumstats::format_sumstats(sumstats_dt,
+                                                    ref_genome = "GRCh37",
+                                                    INFO_filter = 0.9,
+                                                    on_ref_genome = FALSE,
+                                                    strand_ambig_filter = FALSE,
+                                                    bi_allelic_filter = FALSE,
+                                                    allele_flip_check = FALSE,
+                                                    log_folder_ind = FALSE,
+                                                    imputation_ind = FALSE,
+                                                    dbSNP=144,
+                                                    ignore_multi_trait=TRUE,
+                                                    return_data = TRUE)
+      testthat::expect_equal(reformatted$PVAL_ORG,reformatted$P)
+      #test log10
+      sumstats_dt <- data.table::fread(system.file("extdata", "eduAttainOkbay.txt",
+                                                   package = "MungeSumstats"), 
+                                       nThread = 1)
+      sumstats_dt[,LP:=log(Pval,base=10)]
+      data.table::setnames(sumstats_dt,"Pval","Pval_org")
+      reformatted <- MungeSumstats::format_sumstats(sumstats_dt,
+                                                    ref_genome = "GRCh37",
+                                                    INFO_filter = 0.9,
+                                                    on_ref_genome = FALSE,
+                                                    strand_ambig_filter = FALSE,
+                                                    bi_allelic_filter = FALSE,
+                                                    allele_flip_check = FALSE,
+                                                    log_folder_ind = FALSE,
+                                                    imputation_ind = FALSE,
+                                                    dbSNP=144,
+                                                    ignore_multi_trait=TRUE,
+                                                    return_data = TRUE)
+      testthat::expect_equal(reformatted$PVAL_ORG,reformatted$P)
+    }    
+    else{
+      testthat::expect_equal(is_32bit_windows, TRUE)
+      testthat::expect_equal(is_32bit_windows, TRUE)
+    }
+})