LDpred2.R: if --merge-by-rsid is used, ignore --col-bp <arg> (#193)

* LDpred2.R: if `--merge-by-rsid` option is used, ignore `--col-bp` arg
Fixes #191

* should check for NULL, not NA w. argparser

* --merge-by-rsid should be boolean flag

Author: espenhgn

CHANGELOG.md

@@ -39,6 +39,12 @@ If MD5 sum is not listed for a certain release then it means that the container
 
 * Miscellaneous goes here
 
+## [1.3.6] - 2023-08-17
+
+### Fixed
+
+* Ignore LDpred2 `--col-bp <column>` arg in case `--merge-by-rsid` is used
+
 ## [1.3.5] - 2023-08-17
 
 ### Updated

scripts/pgs/LDpred2/ldpred2.R

@@ -34,7 +34,7 @@ par <- add_argument(par, "--col-chr", help="CHR number column", default="CHR", n
 par <- add_argument(par, "--col-snp-id", help="SNP ID (RSID) column", default="SNP", nargs=1)
 par <- add_argument(par, "--col-A1", help="Effective allele column", default="A1", nargs=1)
 par <- add_argument(par, "--col-A2", help="Noneffective allele column", default="A2", nargs=1)
-par <- add_argument(par, "--col-bp", help="SNP position column", default="BP", nargs=1)
+par <- add_argument(par, "--col-bp", help="SNP position column. Will be ignored if --merge-by-rsid flag is used", default="BP", nargs=1)
 par <- add_argument(par, "--col-stat", help="Effect estimate column", default="BETA", nargs=1)
 par <- add_argument(par, "--col-stat-se", help="Effect estimate standard error column", default="SE", nargs=1)
 par <- add_argument(par, "--col-pvalue", help="P-value column", default="P", nargs=1)
@@ -52,7 +52,7 @@ par <- add_argument(par, "--hyper-p-max", help="Maximum (<1) of hyperparameter p
 par <- add_argument(par, "--ldpred-mode", help='Ether "auto" or "inf" (infinitesimal)', default="inf")
 par <- add_argument(par, "--cores", help="Number of CPU cores to use, otherwise use the available number of cores minus 1", default=nb_cores())
 par <- add_argument(par, '--set-seed', help="Set a seed for reproducibility", nargs=1)
-par <- add_argument(par, "--merge-by-rsid", help="Merge using rsid (the default is to merge by chr:bp:a1:a2 codes)", nargs=0)
+par <- add_argument(par, "--merge-by-rsid", help="Merge using rsid (the default is to merge by chr:bp:a1:a2 codes).", flag=TRUE)
 
 parsed <- parse_args(par)
 
@@ -86,7 +86,14 @@ colStat <- parsed$col_stat
 colStatSE <- parsed$col_stat_se
 colPValue <- parsed$col_pvalue
 colN <- parsed$col_n
-mergeByRsid <- !is.na(parsed$merge_by_rsid)
+mergeByRsid <- parsed$merge_by_rsid
+
+# unset colBP in case of merging by rsid
+if (mergeByRsid) {
+  cat(paste('The --merge-by-rsid flag is used; --col-bp', colBP, 'will be ignored\n'))
+  colBP <- NULL
+}
+
 # Polygenic score
 nameScore <- parsed$name_score
 # If the PGS is assigned a name, the diagnostic plot (auto mode only)
@@ -113,8 +120,13 @@ setSeed <- parsed$set_seed
 
 # These vectors are used to convert headers in the sumstat files to those
 # used by bigsnpr
-colSumstatsOld <- c(  colChr, colSNPID, colBP, colA1, colA2, colStat, colStatSE)
-colSumstatToGeno <- c("chr",  "rsid",  "pos",  "a1",  "a0",  "beta",  "beta_se")
+if (mergeByRsid) {
+  colSumstatsOld <- c(colChr, colSNPID, colA1, colA2, colStat, colStatSE)
+  colSumstatToGeno <- c("chr",  "rsid",  "a1",  "a0",  "beta",  "beta_se")
+} else {
+  colSumstatsOld <- c(colChr, colSNPID, colBP, colA1, colA2, colStat, colStatSE)
+  colSumstatToGeno <- c("chr",  "rsid",  "pos",  "a1",  "a0",  "beta",  "beta_se")
+}
 
 # If the user has requested to merge scores to an existing output file
 if (fileOutputMerge) {

tests/test_LDpred2/scripts/extended.sh

@@ -1,7 +1,7 @@
 LDP="$RSCRIPT $DIR_SCRIPTS/ldpred2.R --file-keep-snps $fileKeepSNPS \
   --merge-by-rsid \
   --col-stat beta --col-stat-se beta_se \
-  --col-snp-id rsid --col-bp pos --col-A1 a1 --col-A2 a0 \
+  --col-snp-id rsid --col-bp ignored --col-A1 a1 --col-A2 a0 \
   --out ${fileOut}_imputed.inf"
 
 # Create sumstats with characters in chromosome column
@@ -36,14 +36,14 @@ if [ $? -eq 1 ]; then echo "$dump"; exit; fi
 LDP="$RSCRIPT $DIR_SCRIPTS/ldpred2.R --file-keep-snps $fileKeepSNPS \
   --merge-by-rsid \
   --col-stat beta --col-stat-se beta_se \
-  --col-snp-id rsid --col-chr Chr --col-bp pos --col-A1 a0 --col-A2 a1 \
+  --col-snp-id rsid --col-chr Chr --col-bp ignored --col-A1 a0 --col-A2 a1 \
   --geno-file-rds $fileOutputSNPR --sumstats $fileInputSumStats"
 
-echo "TEST error: Bad ldpred mode"
+echo "Test error: Bad ldpred mode"
 dump=$( { $LDP --ldpred-mode infinite --out $fileOut.inf; } 2>&1 )
 if [ $? -eq 0 ]; then echo "No error received"; exit; fi
 
-echo "TEST error: Bad imputation mode"
+echo "Test error: Bad imputation mode"
 dump=$( { $LDP --ldpred-mode inf --geno-impute bad-mode --out $fileOut.inf; } 2>&1 )
 if [ $? -eq 0 ]; then echo "No error received"; exit; fi
 

version/version.py

@@ -2,7 +2,7 @@
 _MINOR = "3"
 # On main and in a nightly release the patch should be one ahead of the last
 # released build.
-_PATCH = "5"
+_PATCH = "6"
 # This is mainly for nightly builds which have the suffix ".dev$DATE". See
 # https://semver.org/#is-v123-a-semantic-version for the semantics.
 _SUFFIX = ""