.gitignore

+
+# /
+/.project
+/autom4te.cache
+/aclocal.m4
+/Makefile.in
+/configure
+
+# /checks/
+/checks/Makefile
+/checks/Makefile.in
+
+# /checks/R-tests/
+/checks/R-tests/Makefile.in
+
+# /doc/
+/doc/Makefile.in
+/doc/doxygen-docs
+
+# /examples/
+/examples/Makefile.in
+
+# /src/
+/src/Makefile.in
+/src/config.h.in
+.Rhistory

.travis.yml

+dist: trusty
+language: cpp
+script:
+  - autoreconf -i
+  - ./configure --enable-dev-options --disable-silent-rules
+  - make
+  - make check
+compiler:
+  - clang
+  - gcc
+before_install:
+  - sudo apt-get update -qq
+  - sudo apt-get install -qq libeigen3-dev

README.org

+* ProbABEL
+  [[https://api.travis-ci.org/GenABEL-Project/ProbABEL.svg?branch=master]]
+
+  ProbABEL is a tool for genome-wide association analysis of imputed
+  genetic data. It is part of the GenABEL Suite of tools developed by
+  within [[http://www.genabel.org][the GenABEL Project]].
+** Installation
+   Please see the =INSTALL= file in the =doc= directory for more
+   information on how to install ProbABEL.
+** Support
+   For end-user support to ProbABEL please have a look at our forum at
+   http://forum.genabel.org. If you'd like to be kept up-to-date with
+   the latest releases of the tools in the GenABEL suite, please
+   subscribe to the =genabel-announce= mailing list [[https://r-forge.r-project.org/mail/?group_id=505][here]].
+
+   (Potential) Developers are encouraged to sign up for [[https://r-forge.r-project.org/mail/?group_id=505][the
+   =genabel-devel= mailing list on R-forge]].
+** Licence
+   ProbABEL is licensed under the GNU Public Licence v2 or later.
+
+   #+begin_example
+   This program is free software; you can redistribute it and/or
+   modify it under the terms of the GNU General Public License
+   as published by the Free Software Foundation; either version 2
+   of the License, or (at your option) any later version.
+
+   This program is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+   GNU General Public License for more details.
+
+   You should have received a copy of the GNU General Public License
+   along with this program; if not, write to the Free Software
+   Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
+   #+end_example

checks/Makefile.am

@@ -44,19 +44,34 @@ check_SCRIPTS =
 
 if BUILD_palinear
 check_SCRIPTS += check_probabel_chunk.sh check_dose_input.sh
-check_SCRIPTS += test_qt.sh test_mms.sh
+check_SCRIPTS += test_qt.sh test_mms.sh test_flipmaf_palinear.sh
 endif
 if BUILD_palogist
-check_SCRIPTS += test_bt.sh
+check_SCRIPTS += test_bt.sh test_flipmaf_palogist.sh
 endif
 if BUILD_pacoxph
-check_SCRIPTS += test_cox.sh test_cox_nocovar.sh
+check_SCRIPTS += test_cox.sh			\
+ test_cox_nocovar.sh				\
+ test_cox_interactions.sh			\
+ test_flipmaf_pacoxph.sh
 endif
 
 AM_TESTS_ENVIRONMENT= \
+	AWK=$(AWK) \
+	SED=$(SED) \
 	PA_BINDIR=$(top_builddir)/src
 TESTS = $(check_SCRIPTS)
 
+## Dependencies for the test scripts
+TEST_EXTENSIONS=.sh
+test_flipmaf_palinear.log: test_qt.log
+test_flipmaf_palogist.log: test_bt.log
+test_flipmaf_pacoxph.log: test_cox.log
+
+## The pacoxph nocovar test doesn't run correctly until bug #1266 is
+## fixed.
+XFAIL_TESTS = check_pacoxph_nocovar.sh
+
 EXTRA_DIST = $(verified_results) $(input_files) $(check_SCRIPTS)
 
 
@@ -71,6 +86,19 @@ logist_prob_fv_add.out.txt logist_prob_over_domin.out.txt	\
 logist_prob_domin.out.txt logist_prob_fv_recess.out.txt		\
 logist_prob_fv_over_domin.out.txt logist_prob_recess.out.txt
 
+cleanfiles_bt_flipmaf = logist_flipmaf_add.out.txt	\
+logist_flipmaf_fv_add.out.txt				\
+logist_prob_flipmaf_2df.out.txt				\
+logist_prob_flipmaf_fv_2df.out.txt			\
+logist_prob_flipmaf_add.out.txt				\
+logist_prob_flipmaf_fv_add.out.txt			\
+logist_prob_flipmaf_domin.out.txt			\
+logist_prob_flipmaf_fv_domin.out.txt			\
+logist_prob_flipmaf_over_domin.out.txt			\
+logist_prob_flipmaf_fv_over_domin.out.txt		\
+logist_prob_flipmaf_recess.out.txt			\
+logist_prob_flipmaf_fv_recess.out.txt
+
 cleanfiles_qt = linear_base_add.out.txt linear_base_fv_add.out.txt	\
 linear_allcov_add.out.txt linear_allcov_fv_add.out.txt			\
 linear_int1_add.out.txt linear_int1_fv_add.out.txt			\
@@ -112,6 +140,19 @@ linear_ngp2_robust_int1_fv_2df.out.txt					\
 linear_ngp2_allcov_fv_2df.out.txt linear_ngp2_robust_fv_add.out.txt	\
 linear_ngp2_allcov_2df.out.txt
 
+cleanfiles_qt_flipmaf = linear_flipmaf_add.out.txt	\
+linear_flipmaf_fv_add.out.txt				\
+linear_ngp2_flipmaf_2df.out.txt				\
+linear_ngp2_flipmaf_add.out.txt				\
+linear_ngp2_flipmaf_domin.out.txt			\
+linear_ngp2_flipmaf_over_domin.out.txt			\
+linear_ngp2_flipmaf_recess.out.txt			\
+linear_ngp2_flipmaf_fv_2df.out.txt			\
+linear_ngp2_flipmaf_fv_add.out.txt			\
+linear_ngp2_flipmaf_fv_domin.out.txt			\
+linear_ngp2_flipmaf_fv_over_domin.out.txt		\
+linear_ngp2_flipmaf_fv_recess.out.txt
+
 cleanfiles_mms = mmscore_dose_add.out.txt mmscore_dose_fv_add.out.txt	\
 mmscore_prob_fv_over_domin.out.txt mmscore_prob_fv_domin.out.txt	\
 mmscore_prob_over_domin.out.txt mmscore_prob_fv_add.out.txt		\
@@ -119,12 +160,24 @@ mmscore_prob_fv_recess.out.txt mmscore_prob_domin.out.txt		\
 mmscore_prob_recess.out.txt mmscore_prob_add.out.txt			\
 mmscore_prob_2df.out.txt mmscore_prob_fv_2df.out.txt
 
-cleanfiles_cox = coxph_dose_add.out.txt coxph_dose_fv_add.out.txt	\
-coxph_prob_fv_over_domin.out.txt coxph_prob_fv_domin.out.txt		\
-coxph_prob_over_domin.out.txt coxph_prob_fv_add.out.txt			\
-coxph_prob_fv_recess.out.txt coxph_prob_domin.out.txt			\
-coxph_prob_recess.out.txt coxph_prob_add.out.txt			\
-coxph_prob_2df.out.txt coxph_prob_fv_2df.out.txt
+cleanfiles_cox = pacoxph_nocovar_add.out.txt			\
+coxph_data.txt							\
+coxph_dose_add.out.txt coxph_dose_fv_add.out.txt		\
+coxph_prob_fv_over_domin.out.txt coxph_prob_fv_domin.out.txt	\
+coxph_prob_over_domin.out.txt coxph_prob_fv_add.out.txt		\
+coxph_prob_fv_recess.out.txt coxph_prob_domin.out.txt		\
+coxph_prob_recess.out.txt coxph_prob_add.out.txt		\
+coxph_prob_2df.out.txt coxph_prob_fv_2df.out.txt		\
+coxph_prob_allcov_recess.out.txt				\
+coxph_prob_allcov_fv_2df.out.txt				\
+coxph_prob_allcov_add.out.txt					\
+coxph_prob_flipmaf_allcov_over_domin.out.txt			\
+coxph_prob_allcov_fv_over_domin.out.txt				\
+coxph_prob_flipmaf_allcov_fv_domin.out.txt			\
+coxph_prob_flipmaf_allcov_fv_add.out.txt			\
+coxph_prob_flipmaf_allcov_add.out.txt				\
+coxph_prob_allcov_fv_domin.out.txt				\
+coxph_prob_allcov_fv_recess.out.txt
 
 cleanfiles_cox_nocovar = coxph_dose_nocovar_add.out.txt	\
 coxph_dose_nocovar_fv_add.out.txt			\
@@ -139,6 +192,50 @@ coxph_prob_nocovar_add.out.txt				\
 coxph_prob_nocovar_2df.out.txt				\
 coxph_prob_nocovar_fv_2df.out.txt
 
+cleanfiles_cox_interaction = coxph_dose_int1_add.out.txt	\
+coxph_prob_int1_over_domin.out.txt				\
+coxph_prob_int1_domin.out.txt					\
+coxph_prob_int1_recess.out.txt					\
+coxph_prob_int1_add.out.txt					\
+coxph_prob_int1_2df.out.txt					\
+coxph_dose_int2_add.out.txt					\
+coxph_prob_int2_over_domin.out.txt				\
+coxph_prob_int2_domin.out.txt					\
+coxph_prob_int2_recess.out.txt					\
+coxph_prob_int2_add.out.txt					\
+coxph_prob_int2_2df.out.txt					\
+coxph_dose_int3_add.out.txt					\
+coxph_prob_int3_over_domin.out.txt				\
+coxph_prob_int3_domin.out.txt					\
+coxph_prob_int3_recess.out.txt					\
+coxph_prob_int3_add.out.txt					\
+coxph_prob_int3_2df.out.txt
+
+
+cleanfiles_cox_flipmaf = coxph_flipmaf_add.out.txt	\
+coxph_flipmaf_fv_add.out.txt				\
+coxph_prob_flipmaf_2df.out.txt				\
+coxph_prob_flipmaf_add.out.txt				\
+coxph_prob_flipmaf_domin.out.txt			\
+coxph_prob_flipmaf_over_domin.out.txt			\
+coxph_prob_flipmaf_recess.out.txt			\
+coxph_prob_flipmaf_fv_2df.out.txt			\
+coxph_prob_flipmaf_fv_add.out.txt			\
+coxph_prob_flipmaf_fv_domin.out.txt			\
+coxph_prob_flipmaf_fv_over_domin.out.txt		\
+coxph_prob_flipmaf_fv_recess.out.txt			\
+coxph_prob_allcov_over_domin.out.txt			\
+coxph_prob_allcov_fv_add.out.txt			\
+coxph_prob_allcov_2df.out.txt				\
+coxph_prob_flipmaf_allcov_fv_over_domin.out.txt		\
+coxph_prob_flipmaf_allcov_2df.out.txt			\
+coxph_prob_flipmaf_allcov_fv_2df.out.txt		\
+coxph_prob_flipmaf_allcov_fv_recess.out.txt		\
+coxph_prob_allcov_domin.out.txt				\
+coxph_prob_flipmaf_allcov_recess.out.txt		\
+coxph_prob_flipmaf_allcov_domin.out.txt
+
+
 dose_files = chr1.chunk1.dose chr1.chunk2.dose chr2.chunk1.dose	\
  chr2.chunk2.dose chr1.dose chr2.dose
 
@@ -159,6 +256,8 @@ cleanfiles_probabel_check = height.PHE $(dose_files) $(prob_files)	\
 
 cleanfiles_dose_check = test.mldose
 
-CLEANFILES =  $(cleanfiles_probabel_check) $(cleanfiles_dose_check)	\
- $(cleanfiles_bt) $(cleanfiles_qt) $(cleanfiles_mms) $(cleanfiles_cox)  \
- $(cleanfiles_cox_nocovar)
+CLEANFILES = $(cleanfiles_probabel_check) $(cleanfiles_dose_check)	\
+ $(cleanfiles_bt) $(cleanfiles_bt_flipmaf) $(cleanfiles_qt)		\
+ $(cleanfiles_qt_flipmaf) $(cleanfiles_mms) $(cleanfiles_cox)		\
+ $(cleanfiles_cox_nocovar) $(cleanfiles_cox_interaction)		\
+ $(cleanfiles_cox_flipmaf)

checks/R-tests/Makefile.am

@@ -29,11 +29,9 @@ TESTS = $(check_SCRIPTS)
 
 ## The palogist R test still doesn't run correctly.
 XFAIL_TESTS = run_R_test_palogist.sh
-
-## The pacoxph R test fails on SNP 6 when EIGEN is not enabled.
-if !WITH_EIGEN
+# The pacoxph R test started to fail after the R survival package was
+# upgraded to v2.39.2 (see issue #27)
 XFAIL_TESTS += run_R_test_pacox.sh
-endif
 
 EXTRA_DIST = $(check_SCRIPTS) $(R_test_files)
 
@@ -102,7 +100,17 @@ coxph_prob_fv_over_domin.out.txt coxph_prob_fv_domin.out.txt		\
 coxph_prob_over_domin.out.txt coxph_prob_fv_add.out.txt			\
 coxph_prob_fv_recess.out.txt coxph_prob_domin.out.txt			\
 coxph_prob_recess.out.txt coxph_prob_add.out.txt			\
-coxph_prob_2df.out.txt coxph_prob_fv_2df.out.txt
+coxph_prob_2df.out.txt coxph_prob_fv_2df.out.txt			\
+coxph_prob_allcov_recess.out.txt					\
+coxph_prob_allcov_fv_2df.out.txt					\
+coxph_prob_allcov_2df.out.txt						\
+coxph_prob_allcov_domin.out.txt						\
+coxph_prob_allcov_fv_over_domin.out.txt					\
+coxph_prob_allcov_over_domin.out.txt					\
+coxph_prob_allcov_fv_add.out.txt					\
+coxph_prob_allcov_fv_recess.out.txt					\
+coxph_prob_allcov_add.out.txt						\
+coxph_prob_allcov_fv_domin.out.txt
 
 cleanfiles_cox_nocovar = coxph_dose_nocovar_add.out.txt	\
 coxph_dose_nocovar_fv_add.out.txt			\

checks/R-tests/initial_checks.R

@@ -11,6 +11,13 @@ tol <- 1e-5
 ## look for the variables freq and poly.
 rsq.thresh <- 1e-16
 
+## A function that prints text with empty space following it. Used to
+## print what is tested (to be followed with an OK if all is OK).
+prnt <- function(string) {
+    cat(sprintf("%-50s", string))
+}
+
+
 ####
 ## load the data
 ####
@@ -28,7 +35,7 @@ dose <- read.table(paste0(inputfiles.path, "test.mldose"),
 idNames   <- dose[, 1]
 idNames   <- sub("[0-9]+->", "", idNames)
 dose[, 1] <- idNames
-cat("Dose: check consistency of names\t\t")
+prnt("Dose: check consistency of names")
 stopifnot( all.equal(dose[, 1], pheno[, 1], tol) )
 cat("OK\n")
 
@@ -39,7 +46,7 @@ prob <- read.table(paste0(inputfiles.path, "test.mlprob"),
 idNames   <- prob[, 1]
 idNames   <- sub("[0-9]+->", "", idNames)
 prob[, 1] <- idNames
-cat("Prob: check consistency of names\t\t")
+prnt("Prob: check consistency of names")
 stopifnot( all.equal(prob[, 1], pheno[, 1], tol) )
 cat("OK\n")
 
@@ -50,7 +57,7 @@ for (i in 3:dim(dose)[2]) {
         indexHet <- indexHom + 1
         doseFromProb[, i] <- prob[, indexHom] * 2 + prob[, indexHet]
 }
-cat("Check consistency dose <-> prob gtdata\t\t")
+prnt("Check consistency dose <-> prob gtdata")
 stopifnot( all.equal(dose[, 3:ncol(dose)],
                      as.data.frame(doseFromProb)[,3:ncol(doseFromProb)],
                      tol=tol )
@@ -65,26 +72,22 @@ Rsq <- read.table(paste0(inputfiles.path, "test.mlinfo"),
 
 
 ## Define column names of the various ProbABEL output file headers
-colsAddDose <- c("Rsq",
-                 "beta_SNP_add",
-                 "sebeta_SNP_add",
-                 "chi2_SNP")
-colsAddProb <- c("Rsq",
-                 "beta_SNP_addA1",
-                 "sebeta_SNP_addA1",
-                 "chi2_SNP_A1")
+colsAdd <- c("Rsq",
+             "beta_SNP_addA1",
+             "sebeta_SNP_addA1",
+             "chi2_SNP_add")
 colsDom <- c("Rsq",
              "beta_SNP_domA1",
              "sebeta_SNP_domA1",
-             "chi2_SNP_domA1")
+             "chi2_SNP_dom")
 colsRec <- c("Rsq",
              "beta_SNP_recA1",
              "sebeta_SNP_recA1",
-             "chi2_SNP_recA1")
+             "chi2_SNP_rec")
 colsOdom <-c("Rsq",
              "beta_SNP_odomA1",
              "sebeta_SNP_odomA1",
-             "chi2_SNP_odomA1")
+             "chi2_SNP_odom")
 cols2df <- c("Rsq",
              "beta_SNP_A1A2",
              "sebeta_SNP_A1A2",

checks/R-tests/run_model_coxph.R

@@ -10,10 +10,12 @@
 ##' @param snpcomponent2 String telling how the second SNP term is
 ##' defined (only used in the 2 df model). By default this term is
 ##' constant ("1")
+##' @param verbose Logical indicating whether extra debug messages
+##'     should be printed. Default: FALSE
 ##' @return A data frame containing the coefficients from the
 ##' regression analysis and some other variables, such that this
 ##' output can be compared to the ProbABEL output.
-##' @author L.C. Larsen
+##' @author L.C. Karssen
 run.model <- function(model0.txt, model.txt,
                       snpcomponent1, snpcomponent2="1",
                       verbose=FALSE) {
@@ -29,6 +31,7 @@ run.model <- function(model0.txt, model.txt,
 
     for (i in 3:dim(dose)[2]) {
         if (verbose) print(paste("------- new iteration: i =", i))
+
         indexHom <- 3 + ( i - 3 ) * 2
         indexHet <- indexHom + 1
         snp1     <- eval(parse(text=snpcomponent1))
@@ -37,6 +40,11 @@ run.model <- function(model0.txt, model.txt,
         mu       <- rep(1., length(snp)) # Add a constant mean
 
         noNA    <- which( !is.na(snp) )
+        if (verbose) {
+            print(paste("Null model:", model0.txt))
+            print(paste("Alt. model:", model.txt))
+        }
+
         model.0 <- eval(parse(text=model0.txt))
 
         ## Check the imputation R^2, if below threshold ProbABEL will
@@ -75,12 +83,12 @@ run.model <- function(model0.txt, model.txt,
 
         if (verbose) print(model)
 
-        if ( grepl("Inf", model[[2]]$message) |
+        if (grepl("Inf", model[[2]]$message) |
             grepl("infinite", model[[2]]$message) |
             grepl("iterations", model[[2]]$message) |
-            ( grepl("singular", model[[2]]$message) &
-                 (regexpr("variable 1$", model[[2]]$message)[[1]] != 33)
-             )
+            (grepl("singular", model[[2]]$message) &
+             (regexpr("variable 1$", model[[2]]$message)[[1]] != 33)
+            )
             ) {
             ## The model did not converge or some other
             ## errors/warnings occurred, fill the coefficients with

checks/R-tests/run_models_in_R_pacox.R

@@ -25,15 +25,15 @@ source(paste0(srcdir, "initial_checks.R"))
 ####
 ## Run ProbABEL to get the output data we want to compare/verify
 ####
-cat("Running ProbABEL...\t\t\t\t")
+prnt("Running ProbABEL...")
 tmp <- system(paste0("bash ", tests.path, "test_cox.sh 2> /dev/null"),
               intern=TRUE)
 cat("OK\n")
 
 dose.add.PA <- read.table("coxph_dose_add.out.txt",
-                          head=TRUE)[, colsAddDose]
+                          head=TRUE)[, colsAdd]
 prob.add.PA <- read.table("coxph_prob_add.out.txt",
-                          head=TRUE)[, colsAddProb]
+                          head=TRUE)[, colsAdd]
 prob.dom.PA <- read.table("coxph_prob_domin.out.txt",
                           head=TRUE)[, colsDom]
 prob.rec.PA <- read.table("coxph_prob_recess.out.txt",
@@ -48,7 +48,7 @@ prob.2df.PA <- read.table("coxph_prob_2df.out.txt",
 ####
 attach(pheno)
 
-cat("Comparing R output with ProbABEL output\t\t")
+cat("Comparing R output with ProbABEL output:\n")
 
 source(paste0(srcdir, "run_model_coxph.R"))
 
@@ -57,48 +57,52 @@ model.fn.0 <-
 model.fn <- "coxph( Surv(fupt_chd, chd) ~ mu + sex + age + othercov + snp1 )"
 
 ## Additive model, dosages
+prnt(" additive (dosages)")
 snpdose <- "dose[, i]"
 dose.add.R <- run.model(model.fn.0, model.fn, snpdose)
-colnames(dose.add.R) <- colsAddDose
+colnames(dose.add.R) <- colsAdd
 rownames(dose.add.R) <- NULL
 stopifnot( all.equal(dose.add.PA, dose.add.R, tol=tol) )
-cat("additive ")
-
+cat("OK\n")
 
 ## Additive model, probabilities
+prnt(" additive (probabilities)")
 snpprob <- "doseFromProb[, i]"
 prob.add.R <- run.model(model.fn.0, model.fn, snpprob)
-colnames(prob.add.R) <- colsAddProb
+colnames(prob.add.R) <- colsAdd
 rownames(prob.add.R) <- NULL
 stopifnot( all.equal(prob.add.PA, prob.add.R, tol=tol) )
-cat("additive ")
+cat("OK\n")
 
 ## dominant model
+prnt(" dominant")
 snpprob <- "prob[, indexHom] + prob[, indexHet]"
 prob.dom.R <- run.model(model.fn.0, model.fn, snpprob)
 colnames(prob.dom.R) <- colsDom
 rownames(prob.dom.R) <- NULL
 stopifnot( all.equal(prob.dom.PA, prob.dom.R, tol=tol) )
-cat("dominant ")
+cat("OK\n")
 
 ## recessive model
+prnt(" recessive")
 snpprob <- "prob[, indexHom]"
 prob.rec.R <- run.model(model.fn.0, model.fn, snpprob)
 colnames(prob.rec.R) <- colsRec
 rownames(prob.rec.R) <- NULL
 stopifnot( all.equal(prob.rec.PA, prob.rec.R, tol=tol) )
-cat("recessive ")
+cat("OK\n")
 
 ## over-dominant model
+prnt(" overdominant")
 snpprob <- "prob[, indexHet]"
 prob.odom.R <- run.model(model.fn.0, model.fn, snpprob)
 colnames(prob.odom.R) <- colsOdom
 rownames(prob.odom.R) <- NULL
 stopifnot( all.equal(prob.odom.PA, prob.odom.R, tol=tol) )
-cat("overdominant ")
-
+cat("OK\n")
 
 ## 2df model
+prnt(" 2df")
 model.fn <-
     "coxph( Surv(fupt_chd, chd) ~ sex + age + othercov + snp1 + snp2 )"
 snpd1 <- "prob[, indexHet]"
@@ -107,20 +111,20 @@ prob.2df.R <- run.model(model.fn.0, model.fn, snpd1, snpd2)
 colnames(prob.2df.R) <- cols2df
 rownames(prob.2df.R) <- NULL
 stopifnot( all.equal(prob.2df.PA, prob.2df.R, tol=tol) )
-cat("2df\n")
+cat("OK\n")
 
 
 
 cat("\nRun the same checks, but without any covariates (see bug #1266)\n")
-cat("Running ProbABEL...\t\t\t\t")
+prnt("Running ProbABEL...")
 tmp <- system(paste0("bash ", tests.path, "test_cox_nocovar.sh 2> /dev/null"),
               intern=TRUE)
 cat("OK\n")
 
 dose.add.PA <- read.table("coxph_dose_nocovar_add.out.txt",
-                          head=TRUE)[, colsAddDose]
+                          head=TRUE)[, colsAdd]
 prob.add.PA <- read.table("coxph_prob_nocovar_add.out.txt",
-                          head=TRUE)[, colsAddProb]
+                          head=TRUE)[, colsAdd]
 prob.dom.PA <- read.table("coxph_prob_nocovar_domin.out.txt",
                           head=TRUE)[, colsDom]
 prob.rec.PA <- read.table("coxph_prob_nocovar_recess.out.txt",
@@ -135,51 +139,57 @@ model.fn.0 <-
 model.fn <- "coxph( Surv(fupt_chd, chd) ~  mu + snp1 )"
 
 
-cat("Comparing R output with ProbABEL output\t\t")
+cat("Comparing R output with ProbABEL output:\n")
 
 ## Additive model, dosages
+prnt(" additive (dosages)")
 snpdose <- "dose[, i]"
 dose.add.R <- run.model(model.fn.0, model.fn, snpdose)
-colnames(dose.add.R) <- colsAddDose
+colnames(dose.add.R) <- colsAdd
 rownames(dose.add.R) <- NULL
 stopifnot( all.equal(dose.add.PA, dose.add.R, tol=tol) )
-cat("additive ")
+cat("OK\n")
 
 
 ## Additive model, probabilities
+prnt(" additive (probabilities)")
 snpprob <- "doseFromProb[, i]"
 prob.add.R <- run.model(model.fn.0, model.fn, snpprob)
-colnames(prob.add.R) <- colsAddProb
+colnames(prob.add.R) <- colsAdd
 rownames(prob.add.R) <- NULL
 stopifnot( all.equal(prob.add.PA, prob.add.R, tol=tol) )
-cat("additive ")
+cat("OK\n")
 
 ## dominant model
+prnt(" dominant")
 snpprob <- "prob[, indexHom] + prob[, indexHet]"
 prob.dom.R <- run.model(model.fn.0, model.fn, snpprob)
 colnames(prob.dom.R) <- colsDom
 rownames(prob.dom.R) <- NULL
 stopifnot( all.equal(prob.dom.PA, prob.dom.R, tol=tol) )
-cat("dominant ")
+cat("OK\n")
 
 ## recessive model
+prnt(" recessive")
 snpprob <- "prob[, indexHom]"
 prob.rec.R <- run.model(model.fn.0, model.fn, snpprob)
 colnames(prob.rec.R) <- colsRec
 rownames(prob.rec.R) <- NULL
 stopifnot( all.equal(prob.rec.PA, prob.rec.R, tol=tol) )
-cat("recessive ")
+cat("OK\n")
 
 ## over-dominant model
+prnt(" overdominant")
 snpprob <- "prob[, indexHet]"
 prob.odom.R <- run.model(model.fn.0, model.fn, snpprob)
 colnames(prob.odom.R) <- colsOdom
 rownames(prob.odom.R) <- NULL
 stopifnot( all.equal(prob.odom.PA, prob.odom.R, tol=tol) )
-cat("overdominant ")
+cat("OK\n")
 
 
 ## 2df model
+prnt(" 2df")
 model.fn <-
     "coxph( Surv(fupt_chd, chd) ~ mu+ snp1 + snp2 )"
 snpd1 <- "prob[, indexHet]"
@@ -188,6 +198,4 @@ prob.2df.R <- run.model(model.fn.0, model.fn, snpd1, snpd2)
 colnames(prob.2df.R) <- cols2df
 rownames(prob.2df.R) <- NULL
 stopifnot( all.equal(prob.2df.PA, prob.2df.R, tol=tol) )
-cat("2df\n")
-
-cat("\t\t\t\t\t\tOK\n")
+cat("OK\n")

checks/R-tests/run_models_in_R_palinear.R

@@ -14,15 +14,15 @@ source(paste0(srcdir, "initial_checks.R"))
 ####
 ## Run ProbABEL to get the output data we want to compare/verify
 ####
-cat("Running ProbABEL...\t\t\t\t")
+prnt("Running ProbABEL...")
 tmp <- system(paste0("bash ", tests.path, "test_qt.sh"),
               intern=TRUE)
 cat("OK\n")
 
 dose.add.PA <- read.table("linear_base_add.out.txt",
-                          head=TRUE)[, colsAddDose]
+                          head=TRUE)[, colsAdd]
 prob.add.PA <- read.table("linear_ngp2_add.out.txt",
-                          head=TRUE)[, colsAddProb]
+                          head=TRUE)[, colsAdd]
 prob.dom.PA <- read.table("linear_ngp2_domin.out.txt",
                           head=TRUE)[, colsDom]
 prob.rec.PA <- read.table("linear_ngp2_recess.out.txt",
@@ -36,14 +36,22 @@ prob.2df.PA <- read.table("linear_ngp2_2df.out.txt",
 ## (SNP 6 in the info file). ProbABEL lists them all as 0.0, R lists
 ## them as:
 prob.dom.PA[6, 2:4] <- c(NaN, NaN, 0.0)
-prob.2df.PA[6, 4:5] <- c(NA, NA)
+
+## For the 2df model the last SNP is interchangeable: ProbABEL
+## calculates the beta for the other SNP than R. This causes the beta
+## to have the wrong sign. This part changes the position of the SNP
+## beta (and swaps sign) and SE if beta and the other SE are 0.
+if (sum(abs(prob.2df.PA[6, 2:3])) == 0){
+    prob.2df.PA[6, 2:3] <- c(prob.2df.PA[6, 4] * -1, prob.2df.PA[6, 5])
+    prob.2df.PA[6, 4:5] <- c(NA, NA)
+}
 
 ####
 ## run analysis in R
 ####
 attach(pheno)
 
-cat("Comparing R output with ProbABEL output\t\t")
+cat("Comparing R output with ProbABEL output:\n")
 
 source(paste0(srcdir, "run_model_linear.R"))
 
@@ -51,47 +59,53 @@ model.fn.0 <- "lm( height[noNA] ~ sex[noNA] + age[noNA] )"
 model.fn   <- "lm( height ~ sex + age + snp )"
 
 ## Additive model, dosages
+prnt(" additive (dosages)")
 snpdose <- "dose[, i]"
 dose.add.R <- run.model(model.fn.0, model.fn, snpdose)
-colnames(dose.add.R) <- colsAddDose
+colnames(dose.add.R) <- colsAdd
 rownames(dose.add.R) <- NULL
 stopifnot( all.equal(dose.add.PA, dose.add.R, tol=tol) )
-cat("additive ")
-
+cat("OK\n")
 
 ## Additive model, probabilities
+prnt(" additive (probabilities)")
 snpprob <- "doseFromProb[, i]"
 prob.add.R <- run.model(model.fn.0, model.fn, snpprob)
-colnames(prob.add.R) <- colsAddProb
+colnames(prob.add.R) <- colsAdd
 rownames(prob.add.R) <- NULL
 stopifnot( all.equal(prob.add.PA, prob.add.R, tol=tol) )
-cat("additive ")
+cat("OK\n")
 
 ## dominant model
+prnt(" dominant")
 snpprob <- "prob[, indexHom] + prob[, indexHet]"
 prob.dom.R <- run.model(model.fn.0, model.fn, snpprob)
 colnames(prob.dom.R) <- colsDom
 rownames(prob.dom.R) <- NULL
 stopifnot( all.equal(prob.dom.PA, prob.dom.R, tol=tol) )
-cat("dominant ")
+cat("OK\n")
 
 ## recessive model
+prnt(" recessive")
 snpprob <- "prob[, indexHom]"
 prob.rec.R <- run.model(model.fn.0, model.fn, snpprob)
 colnames(prob.rec.R) <- colsRec
 rownames(prob.rec.R) <- NULL
 stopifnot( all.equal(prob.rec.PA, prob.rec.R, tol=tol) )
-cat("recessive ")
+#print("OK", justify="right")
+cat("OK\n")
 
 ## over-dominant model
+prnt(" overdominant")
 snpprob <- "prob[, indexHet]"
 prob.odom.R <- run.model(model.fn.0, model.fn, snpprob)
 colnames(prob.odom.R) <- colsOdom
 rownames(prob.odom.R) <- NULL
 stopifnot( all.equal(prob.odom.PA, prob.odom.R, tol=tol) )
-cat("overdominant ")
+cat("OK\n")
 
 ## 2df model
+prnt(" 2df")
 prob.2df.R <- data.frame()
 for (i in 3:dim(dose)[2]) {
         indexHom <- 3 + ( i - 3 ) * 2
@@ -124,6 +138,4 @@ for (i in 3:dim(dose)[2]) {
 colnames(prob.2df.R) <- cols2df
 rownames(prob.2df.R) <- NULL
 stopifnot( all.equal(prob.2df.PA, prob.2df.R, tol=tol) )
-cat("2df\n")
-
-cat("\t\t\t\t\t\tOK\n")
+cat("OK\n")

checks/R-tests/run_models_in_R_palogist.R

@@ -14,15 +14,15 @@ source(paste0(srcdir, "initial_checks.R"))
 ####
 ## Run ProbABEL to get the output data we want to compare/verify
 ####
-cat("Running ProbABEL...\t\t\t\t")
+prnt("Running ProbABEL...")
 tmp <- system(paste0("bash ", tests.path, "test_bt.sh"),
               intern=TRUE)
 cat("OK\n")
 
 dose.add.PA <- read.table("logist_add.out.txt",
-                          head=TRUE)[, colsAddDose]
+                          head=TRUE)[, colsAdd]
 prob.add.PA <- read.table("logist_prob_add.out.txt",
-                          head=TRUE)[, colsAddProb]
+                          head=TRUE)[, colsAdd]
 prob.dom.PA <- read.table("logist_prob_domin.out.txt",
                           head=TRUE)[, colsDom]
 prob.rec.PA <- read.table("logist_prob_recess.out.txt",
@@ -41,7 +41,7 @@ prob.dom.PA[6, 4] <- 0.0
 ####
 attach(pheno)
 
-cat("Comparing R output with ProbABEL output\t\t")
+cat("Comparing R output with ProbABEL output\n")
 
 source(paste0(srcdir, "run_model_logist.R"))
 
@@ -50,47 +50,52 @@ model.fn.0 <-
 model.fn  <- "glm( chd ~ sex + age + othercov + snp, family=binomial )"
 
 ## Additive model, dosages
+prnt(" additive (dosages)")
 snpdose <- "dose[, i]"
 dose.add.R <- run.model(model.fn.0, model.fn, snpdose)
-colnames(dose.add.R) <- colsAddDose
+colnames(dose.add.R) <- colsAdd
 rownames(dose.add.R) <- NULL
 stopifnot( all.equal(dose.add.PA, dose.add.R, tol=tol) )
-cat("additive ")
+cat("OK\n")
 
 ## Additive model, probabilities
+prnt(" additive (probabilities)")
 snpprob <- "doseFromProb[, i]"
 prob.add.R <- run.model(model.fn.0, model.fn, snpprob)
-colnames(prob.add.R) <- colsAddProb
+colnames(prob.add.R) <- colsAdd
 rownames(prob.add.R) <- NULL
 stopifnot( all.equal(prob.add.PA, prob.add.R, tol=tol) )
-cat("additive ")
+cat("OK\n")
 
 ## dominant model
+prnt(" dominant")
 snpprob <- "prob[, indexHom] + prob[, indexHet]"
 prob.dom.R <- run.model(model.fn.0, model.fn, snpprob)
 colnames(prob.dom.R) <- colsDom
 rownames(prob.dom.R) <- NULL
 stopifnot( all.equal(prob.dom.PA, prob.dom.R, tol=tol) )
-cat("dominant ")
+cat("OK\n")
 
 ## recessive model
+prnt(" recessive")
 snpprob <- "prob[, indexHom]"
 prob.rec.R <- run.model(model.fn.0, model.fn, snpprob)
 colnames(prob.rec.R) <- colsRec
 rownames(prob.rec.R) <- NULL
 stopifnot( all.equal(prob.rec.PA, prob.rec.R, tol=tol) )
-cat("recessive ")
+cat("OK\n")
 
 ## over-dominant model
+prnt(" overdominant")
 snpprob <- "prob[, indexHet]"
 prob.odom.R <- run.model(model.fn.0, model.fn, snpprob)
 colnames(prob.odom.R) <- colsOdom
 rownames(prob.odom.R) <- NULL
 stopifnot( all.equal(prob.odom.PA, prob.odom.R, tol=tol) )
-cat("overdominant ")
-
+cat("OK\n")
 
 ## 2df model
+prnt(" 2df")
 prob.2df.R <- data.frame()
 for (i in 3:dim(dose)[2]) {
         indexHom <- 3 + ( i - 3 ) * 2
@@ -124,6 +129,4 @@ for (i in 3:dim(dose)[2]) {
 colnames(prob.2df.R) <- cols2df
 rownames(prob.2df.R) <- NULL
 stopifnot( all.equal(prob.2df.PA, prob.2df.R, tol=tol) )
-cat("2df\n")
-
-cat("\t\t\t\t\t\tOK\n")
+cat("OK\n")

checks/check_dose_input.sh

@@ -8,6 +8,9 @@ echo "Checking palinear with dose data without '->'"
 # Exit with error when one of the steps in the script fails
 set -e
 
+if [ -z ${sed} ]; then
+    SED=sed
+fi
 if [ -z ${srcdir} ]; then
     srcdir="."
 fi
@@ -15,7 +18,7 @@ inputdir="${srcdir}/inputfiles/"
 results="${srcdir}/verified_results/"
 outfile="height_base_add.out.txt"
 
-sed 's/^[[:digit:]]*->//' $inputdir/test.mldose > test.mldose
+${SED} 's/^[[:digit:]]*->//' $inputdir/test.mldose > test.mldose
 
 
 ../src/palinear \

checks/check_probabel_chunk.sh

@@ -5,6 +5,13 @@
 
 echo "Testing probabel..."
 
+if [ -z ${AWK} ]; then
+    AWK=awk
+fi
+if [ -z ${sed} ]; then
+    SED=sed
+fi
+
 # Exit with error when one of the steps in the script fails
 set -e
 
@@ -28,7 +35,7 @@ chunksep="_._chunk_._"
 chrsep="_._chr_._"
 
 # ------ Prepare probabel and the config file ------
-sed 's;"./";"../src/";g' $probabel > probabel
+${SED} 's;"./";"../src/";g' $probabel > probabel
 chmod a+x probabel
 cp $probabelcfg probabel_config.cfg
 chmod +w probabel_config.cfg # Need this for make distcheck to work
@@ -50,50 +57,50 @@ prepare_input ()
         # ------------------ No chunks test -------------------
         # Split the dose, prob and info files up into two chromosomes
         # with some  chunks
-        awk '{print $1,$2,$3,$4}'    $dosefile > chr1.dose
-        awk '{print $1,$2,$5,$6,$7}' $dosefile > chr2.dose
+        ${AWK} '{print $1,$2,$3,$4}'    $dosefile > chr1.dose
+        ${AWK} '{print $1,$2,$5,$6,$7}' $dosefile > chr2.dose
 
-        awk '{print $1,$2,$3,$4,$5,$6}'          $probfile > chr1.prob
-        awk '{print $1,$2,$7,$8,$9,$10,$11,$12}' $probfile > chr2.prob
+        ${AWK} '{print $1,$2,$3,$4,$5,$6}'          $probfile > chr1.prob
+        ${AWK} '{print $1,$2,$7,$8,$9,$10,$11,$12}' $probfile > chr2.prob
 
-        sed -n '1,3p' $infofile >  chr1.info
-        sed -n '1p'   $infofile >  chr2.info
-        sed -n '4,6p' $infofile >> chr2.info
+        ${SED} -n '1,3p' $infofile >  chr1.info
+        ${SED} -n '1p'   $infofile >  chr2.info
+        ${SED} -n '4,6p' $infofile >> chr2.info
 
-        sed -n '1,3p' $mapfile > chr1.map
-        sed -n '1p'   $mapfile >  chr2.map
-        sed -n '4,6p' $mapfile >> chr2.map
+        ${SED} -n '1,3p' $mapfile > chr1.map
+        ${SED} -n '1p'   $mapfile >  chr2.map
+        ${SED} -n '4,6p' $mapfile >> chr2.map
 
         WithOrWithout="without"
     elif [ "$1" = "chunk" ]; then
         # ------------------ Chunks test ----------------------
         # Split the dose and info files up into two chromosomes with
         # some chunks
-        awk '{print $1,$2,$3}'    $dosefile > chr1.chunk1.dose
-        awk '{print $1,$2,$4}'    $dosefile > chr1.chunk2.dose
-        awk '{print $1,$2,$5,$6}' $dosefile > chr2.chunk1.dose
-        awk '{print $1,$2,$7}'    $dosefile > chr2.chunk2.dose
-
-        awk '{print $1,$2,$3,$4}'        $probfile > chr1.chunk1.prob
-        awk '{print $1,$2,$5,$6}'        $probfile > chr1.chunk2.prob
-        awk '{print $1,$2,$7,$8,$9,$10}' $probfile > chr2.chunk1.prob
-        awk '{print $1,$2,$11,$12}'      $probfile > chr2.chunk2.prob
-
-        sed -n '1,2p' $infofile >  chr1.chunk1.info
-        sed -n '1p'   $infofile >  chr1.chunk2.info
-        sed -n '3p'   $infofile >> chr1.chunk2.info
-        sed -n '1p'   $infofile >  chr2.chunk1.info
-        sed -n '4,5p' $infofile >> chr2.chunk1.info
-        sed -n '1p'   $infofile >  chr2.chunk2.info
-        sed -n '6p'   $infofile >> chr2.chunk2.info
-
-        sed -n '1,2p' $mapfile >  chr1.chunk1.map
-        sed -n '1p'   $mapfile >  chr1.chunk2.map
-        sed -n '3p'   $mapfile >> chr1.chunk2.map
-        sed -n '1p'   $mapfile >  chr2.chunk1.map
-        sed -n '4,5p' $mapfile >> chr2.chunk1.map
-        sed -n '1p'   $mapfile >  chr2.chunk2.map
-        sed -n '6p'   $mapfile >> chr2.chunk2.map
+        ${AWK} '{print $1,$2,$3}'    $dosefile > chr1.chunk1.dose
+        ${AWK} '{print $1,$2,$4}'    $dosefile > chr1.chunk2.dose
+        ${AWK} '{print $1,$2,$5,$6}' $dosefile > chr2.chunk1.dose
+        ${AWK} '{print $1,$2,$7}'    $dosefile > chr2.chunk2.dose
+
+        ${AWK} '{print $1,$2,$3,$4}'        $probfile > chr1.chunk1.prob
+        ${AWK} '{print $1,$2,$5,$6}'        $probfile > chr1.chunk2.prob
+        ${AWK} '{print $1,$2,$7,$8,$9,$10}' $probfile > chr2.chunk1.prob
+        ${AWK} '{print $1,$2,$11,$12}'      $probfile > chr2.chunk2.prob
+
+        ${SED} -n '1,2p' $infofile >  chr1.chunk1.info
+        ${SED} -n '1p'   $infofile >  chr1.chunk2.info
+        ${SED} -n '3p'   $infofile >> chr1.chunk2.info
+        ${SED} -n '1p'   $infofile >  chr2.chunk1.info
+        ${SED} -n '4,5p' $infofile >> chr2.chunk1.info
+        ${SED} -n '1p'   $infofile >  chr2.chunk2.info
+        ${SED} -n '6p'   $infofile >> chr2.chunk2.info
+
+        ${SED} -n '1,2p' $mapfile >  chr1.chunk1.map
+        ${SED} -n '1p'   $mapfile >  chr1.chunk2.map
+        ${SED} -n '3p'   $mapfile >> chr1.chunk2.map
+        ${SED} -n '1p'   $mapfile >  chr2.chunk1.map
+        ${SED} -n '4,5p' $mapfile >> chr2.chunk1.map
+        ${SED} -n '1p'   $mapfile >  chr2.chunk2.map
+        ${SED} -n '6p'   $mapfile >> chr2.chunk2.map
 
         WithOrWithout="with"
     else

checks/compare_with_old.sh

+echo height_base_add.out.txt && diff height_base_add.out.txt  test_other/height_base_add.out.txt
+echo height_base_fv_add.out.txt && diff height_base_fv_add.out.txt  test_other/height_base_fv_add.out.txt
+echo height_allcov_add.out.txt && diff height_allcov_add.out.txt  test_other/height_allcov_add.out.txt
+echo height_allcov_fv_add.out.txt && diff height_allcov_fv_add.out.txt  test_other/height_allcov_fv_add.out.txt
+echo height_int1_add.out.txt && diff height_int1_add.out.txt  test_other/height_int1_add.out.txt
+echo height_int1_fv_add.out.txt && diff height_int1_fv_add.out.txt  test_other/height_int1_fv_add.out.txt
+echo height_robust_add.out.txt && diff height_robust_add.out.txt  test_other/height_robust_add.out.txt
+echo height_robust_fv_add.out.txt && diff height_robust_fv_add.out.txt  test_other/height_robust_fv_add.out.txt
+echo height_robust_int1_add.out.txt && diff height_robust_int1_add.out.txt  test_other/height_robust_int1_add.out.txt
+echo height_robust_int1_fv_add.out.txt && diff height_robust_int1_fv_add.out.txt  test_other/height_robust_int1_fv_add.out.txt
+echo height_ngp2_add.out.txt && diff height_ngp2_add.out.txt  test_other/height_ngp2_add.out.txt
+echo height_ngp2_fv_add.out.txt && diff height_ngp2_fv_add.out.txt  test_other/height_ngp2_fv_add.out.txt
+echo height_ngp2_domin.out.txt && diff height_ngp2_domin.out.txt  test_other/height_ngp2_domin.out.txt
+echo height_ngp2_fv_domin.out.txt && diff height_ngp2_fv_domin.out.txt  test_other/height_ngp2_fv_domin.out.txt
+echo height_ngp2_over_domin.out.txt && diff height_ngp2_over_domin.out.txt  test_other/height_ngp2_over_domin.out.txt
+echo height_ngp2_fv_over_domin.out.txt && diff height_ngp2_fv_over_domin.out.txt  test_other/height_ngp2_fv_over_domin.out.txt
+echo height_ngp2_recess.out.txt && diff height_ngp2_recess.out.txt  test_other/height_ngp2_recess.out.txt
+echo height_ngp2_fv_recess.out.txt && diff height_ngp2_fv_recess.out.txt  test_other/height_ngp2_fv_recess.out.txt
+echo height_ngp2_2df.out.txt && diff height_ngp2_2df.out.txt  test_other/height_ngp2_2df.out.txt
+echo height_ngp2_fv_2df.out.txt && diff height_ngp2_fv_2df.out.txt  test_other/height_ngp2_fv_2df.out.txt
+echo height_ngp2_allcov_add.out.txt && diff height_ngp2_allcov_add.out.txt  test_other/height_ngp2_allcov_add.out.txt
+echo height_ngp2_allcov_fv_add.out.txt && diff height_ngp2_allcov_fv_add.out.txt  test_other/height_ngp2_allcov_fv_add.out.txt
+echo height_ngp2_allcov_domin.out.txt && diff height_ngp2_allcov_domin.out.txt  test_other/height_ngp2_allcov_domin.out.txt
+echo height_ngp2_allcov_fv_domin.out.txt && diff height_ngp2_allcov_fv_domin.out.txt  test_other/height_ngp2_allcov_fv_domin.out.txt
+echo height_ngp2_allcov_over_domin.out.txt && diff height_ngp2_allcov_over_domin.out.txt  test_other/height_ngp2_allcov_over_domin.out.txt
+echo height_ngp2_allcov_fv_over_domin.out.txt && diff height_ngp2_allcov_fv_over_domin.out.txt  test_other/height_ngp2_allcov_fv_over_domin.out.txt
+echo height_ngp2_allcov_recess.out.txt && diff height_ngp2_allcov_recess.out.txt  test_other/height_ngp2_allcov_recess.out.txt
+echo height_ngp2_allcov_fv_recess.out.txt && diff height_ngp2_allcov_fv_recess.out.txt  test_other/height_ngp2_allcov_fv_recess.out.txt
+echo height_ngp2_allcov_2df.out.txt && diff height_ngp2_allcov_2df.out.txt  test_other/height_ngp2_allcov_2df.out.txt
+echo height_ngp2_allcov_fv_2df.out.txt && diff height_ngp2_allcov_fv_2df.out.txt  test_other/height_ngp2_allcov_fv_2df.out.txt
+echo height_ngp2_int1_add.out.txt && diff height_ngp2_int1_add.out.txt  test_other/height_ngp2_int1_add.out.txt
+echo height_ngp2_int1_fv_add.out.txt && diff height_ngp2_int1_fv_add.out.txt  test_other/height_ngp2_int1_fv_add.out.txt
+echo height_ngp2_int1_domin.out.txt && diff height_ngp2_int1_domin.out.txt  test_other/height_ngp2_int1_domin.out.txt
+echo height_ngp2_int1_fv_domin.out.txt && diff height_ngp2_int1_fv_domin.out.txt  test_other/height_ngp2_int1_fv_domin.out.txt
+echo height_ngp2_int1_over_domin.out.txt && diff height_ngp2_int1_over_domin.out.txt  test_other/height_ngp2_int1_over_domin.out.txt
+echo height_ngp2_int1_fv_over_domin.out.txt && diff height_ngp2_int1_fv_over_domin.out.txt  test_other/height_ngp2_int1_fv_over_domin.out.txt
+echo height_ngp2_int1_recess.out.txt && diff height_ngp2_int1_recess.out.txt  test_other/height_ngp2_int1_recess.out.txt
+echo height_ngp2_int1_fv_recess.out.txt && diff height_ngp2_int1_fv_recess.out.txt  test_other/height_ngp2_int1_fv_recess.out.txt
+echo height_ngp2_int1_2df.out.txt && diff height_ngp2_int1_2df.out.txt  test_other/height_ngp2_int1_2df.out.txt
+echo height_ngp2_int1_fv_2df.out.txt && diff height_ngp2_int1_fv_2df.out.txt  test_other/height_ngp2_int1_fv_2df.out.txt
+echo height_ngp2_robust_add.out.txt && diff height_ngp2_robust_add.out.txt  test_other/height_ngp2_robust_add.out.txt
+echo height_ngp2_robust_fv_add.out.txt && diff height_ngp2_robust_fv_add.out.txt  test_other/height_ngp2_robust_fv_add.out.txt
+echo height_ngp2_robust_domin.out.txt && diff height_ngp2_robust_domin.out.txt  test_other/height_ngp2_robust_domin.out.txt
+echo height_ngp2_robust_fv_domin.out.txt && diff height_ngp2_robust_fv_domin.out.txt  test_other/height_ngp2_robust_fv_domin.out.txt
+echo height_ngp2_robust_over_domin.out.txt && diff height_ngp2_robust_over_domin.out.txt  test_other/height_ngp2_robust_over_domin.out.txt
+echo height_ngp2_robust_fv_over_domin.out.txt && diff height_ngp2_robust_fv_over_domin.out.txt  test_other/height_ngp2_robust_fv_over_domin.out.txt
+echo height_ngp2_robust_recess.out.txt && diff height_ngp2_robust_recess.out.txt  test_other/height_ngp2_robust_recess.out.txt
+echo height_ngp2_robust_fv_recess.out.txt && diff height_ngp2_robust_fv_recess.out.txt  test_other/height_ngp2_robust_fv_recess.out.txt
+echo height_ngp2_robust_2df.out.txt && diff height_ngp2_robust_2df.out.txt  test_other/height_ngp2_robust_2df.out.txt
+echo height_ngp2_robust_fv_2df.out.txt && diff height_ngp2_robust_fv_2df.out.txt  test_other/height_ngp2_robust_fv_2df.out.txt
+echo height_ngp2_robust_int1_add.out.txt && diff height_ngp2_robust_int1_add.out.txt  test_other/height_ngp2_robust_int1_add.out.txt
+echo height_ngp2_robust_int1_fv_add.out.txt && diff height_ngp2_robust_int1_fv_add.out.txt  test_other/height_ngp2_robust_int1_fv_add.out.txt
+echo height_ngp2_robust_int1_domin.out.txt && diff height_ngp2_robust_int1_domin.out.txt  test_other/height_ngp2_robust_int1_domin.out.txt
+echo height_ngp2_robust_int1_fv_domin.out.txt && diff height_ngp2_robust_int1_fv_domin.out.txt  test_other/height_ngp2_robust_int1_fv_domin.out.txt
+echo height_ngp2_robust_int1_over_domin.out.txt && diff height_ngp2_robust_int1_over_domin.out.txt  test_other/height_ngp2_robust_int1_over_domin.out.txt
+echo height_ngp2_robust_int1_fv_over_domin.out.txt && diff height_ngp2_robust_int1_fv_over_domin.out.txt  test_other/height_ngp2_robust_int1_fv_over_domin.out.txt
+echo height_ngp2_robust_int1_recess.out.txt && diff height_ngp2_robust_int1_recess.out.txt  test_other/height_ngp2_robust_int1_recess.out.txt
+echo height_ngp2_robust_int1_fv_recess.out.txt && diff height_ngp2_robust_int1_fv_recess.out.txt  test_other/height_ngp2_robust_int1_fv_recess.out.txt
+echo height_ngp2_robust_int1_2df.out.txt && diff height_ngp2_robust_int1_2df.out.txt  test_other/height_ngp2_robust_int1_2df.out.txt
+echo height_ngp2_robust_int1_fv_2df.out.txt && diff height_ngp2_robust_int1_fv_2df.out.txt  test_other/height_ngp2_robust_int1_fv_2df.out.txt
+echo mmscore_add.out.txt && diff mmscore_add.out.txt  test_other/mmscore_add.out.txt
+echo mmscore_fv_add.out.txt && diff mmscore_fv_add.out.txt test_other/mmscore_fv_add.out.txt
+echo mmscore_prob_add.out.txt && diff mmscore_prob_add.out.txt test_other/mmscore_prob_add.out.txt 
+echo mmscore_prob_fv_add.out.txt && diff mmscore_prob_fv_add.out.txt test_other/mmscore_prob_fv_add.out.txt

checks/compare_with_old_short.sh

+echo height_base_add.out.txt && diff height_base_add.out.txt  test_other/height_base_add.out.txt |wc -l
+echo height_base_fv_add.out.txt && diff height_base_fv_add.out.txt  test_other/height_base_fv_add.out.txt |wc -l
+echo height_allcov_add.out.txt && diff height_allcov_add.out.txt  test_other/height_allcov_add.out.txt |wc -l
+echo height_allcov_fv_add.out.txt && diff height_allcov_fv_add.out.txt  test_other/height_allcov_fv_add.out.txt |wc -l
+echo height_int1_add.out.txt && diff height_int1_add.out.txt  test_other/height_int1_add.out.txt |wc -l
+echo height_int1_fv_add.out.txt && diff height_int1_fv_add.out.txt  test_other/height_int1_fv_add.out.txt |wc -l
+echo height_robust_add.out.txt && diff height_robust_add.out.txt  test_other/height_robust_add.out.txt |wc -l
+echo height_robust_fv_add.out.txt && diff height_robust_fv_add.out.txt  test_other/height_robust_fv_add.out.txt |wc -l
+echo height_robust_int1_add.out.txt && diff height_robust_int1_add.out.txt  test_other/height_robust_int1_add.out.txt |wc -l
+echo height_robust_int1_fv_add.out.txt && diff height_robust_int1_fv_add.out.txt  test_other/height_robust_int1_fv_add.out.txt |wc -l
+echo height_ngp2_add.out.txt && diff height_ngp2_add.out.txt  test_other/height_ngp2_add.out.txt |wc -l
+echo height_ngp2_fv_add.out.txt && diff height_ngp2_fv_add.out.txt  test_other/height_ngp2_fv_add.out.txt |wc -l
+echo height_ngp2_domin.out.txt && diff height_ngp2_domin.out.txt  test_other/height_ngp2_domin.out.txt |wc -l
+echo height_ngp2_fv_domin.out.txt && diff height_ngp2_fv_domin.out.txt  test_other/height_ngp2_fv_domin.out.txt |wc -l
+echo height_ngp2_over_domin.out.txt && diff height_ngp2_over_domin.out.txt  test_other/height_ngp2_over_domin.out.txt |wc -l
+echo height_ngp2_fv_over_domin.out.txt && diff height_ngp2_fv_over_domin.out.txt  test_other/height_ngp2_fv_over_domin.out.txt |wc -l
+echo height_ngp2_recess.out.txt && diff height_ngp2_recess.out.txt  test_other/height_ngp2_recess.out.txt |wc -l
+echo height_ngp2_fv_recess.out.txt && diff height_ngp2_fv_recess.out.txt  test_other/height_ngp2_fv_recess.out.txt |wc -l
+echo height_ngp2_2df.out.txt && diff height_ngp2_2df.out.txt  test_other/height_ngp2_2df.out.txt |wc -l
+echo height_ngp2_fv_2df.out.txt && diff height_ngp2_fv_2df.out.txt  test_other/height_ngp2_fv_2df.out.txt |wc -l
+echo height_ngp2_allcov_add.out.txt && diff height_ngp2_allcov_add.out.txt  test_other/height_ngp2_allcov_add.out.txt |wc -l
+echo height_ngp2_allcov_fv_add.out.txt && diff height_ngp2_allcov_fv_add.out.txt  test_other/height_ngp2_allcov_fv_add.out.txt |wc -l
+echo height_ngp2_allcov_domin.out.txt && diff height_ngp2_allcov_domin.out.txt  test_other/height_ngp2_allcov_domin.out.txt |wc -l
+echo height_ngp2_allcov_fv_domin.out.txt && diff height_ngp2_allcov_fv_domin.out.txt  test_other/height_ngp2_allcov_fv_domin.out.txt |wc -l
+echo height_ngp2_allcov_over_domin.out.txt && diff height_ngp2_allcov_over_domin.out.txt  test_other/height_ngp2_allcov_over_domin.out.txt |wc -l
+echo height_ngp2_allcov_fv_over_domin.out.txt && diff height_ngp2_allcov_fv_over_domin.out.txt  test_other/height_ngp2_allcov_fv_over_domin.out.txt |wc -l
+echo height_ngp2_allcov_recess.out.txt && diff height_ngp2_allcov_recess.out.txt  test_other/height_ngp2_allcov_recess.out.txt |wc -l
+echo height_ngp2_allcov_fv_recess.out.txt && diff height_ngp2_allcov_fv_recess.out.txt  test_other/height_ngp2_allcov_fv_recess.out.txt |wc -l
+echo height_ngp2_allcov_2df.out.txt && diff height_ngp2_allcov_2df.out.txt  test_other/height_ngp2_allcov_2df.out.txt |wc -l
+echo height_ngp2_allcov_fv_2df.out.txt && diff height_ngp2_allcov_fv_2df.out.txt  test_other/height_ngp2_allcov_fv_2df.out.txt |wc -l
+echo height_ngp2_int1_add.out.txt && diff height_ngp2_int1_add.out.txt  test_other/height_ngp2_int1_add.out.txt |wc -l
+echo height_ngp2_int1_fv_add.out.txt && diff height_ngp2_int1_fv_add.out.txt  test_other/height_ngp2_int1_fv_add.out.txt |wc -l
+echo height_ngp2_int1_domin.out.txt && diff height_ngp2_int1_domin.out.txt  test_other/height_ngp2_int1_domin.out.txt |wc -l
+echo height_ngp2_int1_fv_domin.out.txt && diff height_ngp2_int1_fv_domin.out.txt  test_other/height_ngp2_int1_fv_domin.out.txt |wc -l
+echo height_ngp2_int1_over_domin.out.txt && diff height_ngp2_int1_over_domin.out.txt  test_other/height_ngp2_int1_over_domin.out.txt |wc -l
+echo height_ngp2_int1_fv_over_domin.out.txt && diff height_ngp2_int1_fv_over_domin.out.txt  test_other/height_ngp2_int1_fv_over_domin.out.txt |wc -l
+echo height_ngp2_int1_recess.out.txt && diff height_ngp2_int1_recess.out.txt  test_other/height_ngp2_int1_recess.out.txt |wc -l
+echo height_ngp2_int1_fv_recess.out.txt && diff height_ngp2_int1_fv_recess.out.txt  test_other/height_ngp2_int1_fv_recess.out.txt |wc -l
+echo height_ngp2_int1_2df.out.txt && diff height_ngp2_int1_2df.out.txt  test_other/height_ngp2_int1_2df.out.txt |wc -l
+echo height_ngp2_int1_fv_2df.out.txt && diff height_ngp2_int1_fv_2df.out.txt  test_other/height_ngp2_int1_fv_2df.out.txt |wc -l
+echo height_ngp2_robust_add.out.txt && diff height_ngp2_robust_add.out.txt  test_other/height_ngp2_robust_add.out.txt |wc -l
+echo height_ngp2_robust_fv_add.out.txt && diff height_ngp2_robust_fv_add.out.txt  test_other/height_ngp2_robust_fv_add.out.txt |wc -l
+echo height_ngp2_robust_domin.out.txt && diff height_ngp2_robust_domin.out.txt  test_other/height_ngp2_robust_domin.out.txt |wc -l
+echo height_ngp2_robust_fv_domin.out.txt && diff height_ngp2_robust_fv_domin.out.txt  test_other/height_ngp2_robust_fv_domin.out.txt |wc -l
+echo height_ngp2_robust_over_domin.out.txt && diff height_ngp2_robust_over_domin.out.txt  test_other/height_ngp2_robust_over_domin.out.txt |wc -l
+echo height_ngp2_robust_fv_over_domin.out.txt && diff height_ngp2_robust_fv_over_domin.out.txt  test_other/height_ngp2_robust_fv_over_domin.out.txt |wc -l
+echo height_ngp2_robust_recess.out.txt && diff height_ngp2_robust_recess.out.txt  test_other/height_ngp2_robust_recess.out.txt |wc -l
+echo height_ngp2_robust_fv_recess.out.txt && diff height_ngp2_robust_fv_recess.out.txt  test_other/height_ngp2_robust_fv_recess.out.txt |wc -l
+echo height_ngp2_robust_2df.out.txt && diff height_ngp2_robust_2df.out.txt  test_other/height_ngp2_robust_2df.out.txt |wc -l
+echo height_ngp2_robust_fv_2df.out.txt && diff height_ngp2_robust_fv_2df.out.txt  test_other/height_ngp2_robust_fv_2df.out.txt |wc -l
+echo height_ngp2_robust_int1_add.out.txt && diff height_ngp2_robust_int1_add.out.txt  test_other/height_ngp2_robust_int1_add.out.txt |wc -l
+echo height_ngp2_robust_int1_fv_add.out.txt && diff height_ngp2_robust_int1_fv_add.out.txt  test_other/height_ngp2_robust_int1_fv_add.out.txt |wc -l
+echo height_ngp2_robust_int1_domin.out.txt && diff height_ngp2_robust_int1_domin.out.txt  test_other/height_ngp2_robust_int1_domin.out.txt |wc -l
+echo height_ngp2_robust_int1_fv_domin.out.txt && diff height_ngp2_robust_int1_fv_domin.out.txt  test_other/height_ngp2_robust_int1_fv_domin.out.txt |wc -l
+echo height_ngp2_robust_int1_over_domin.out.txt && diff height_ngp2_robust_int1_over_domin.out.txt  test_other/height_ngp2_robust_int1_over_domin.out.txt |wc -l
+echo height_ngp2_robust_int1_fv_over_domin.out.txt && diff height_ngp2_robust_int1_fv_over_domin.out.txt  test_other/height_ngp2_robust_int1_fv_over_domin.out.txt |wc -l
+echo height_ngp2_robust_int1_recess.out.txt && diff height_ngp2_robust_int1_recess.out.txt  test_other/height_ngp2_robust_int1_recess.out.txt |wc -l
+echo height_ngp2_robust_int1_fv_recess.out.txt && diff height_ngp2_robust_int1_fv_recess.out.txt  test_other/height_ngp2_robust_int1_fv_recess.out.txt |wc -l
+echo height_ngp2_robust_int1_2df.out.txt && diff height_ngp2_robust_int1_2df.out.txt  test_other/height_ngp2_robust_int1_2df.out.txt |wc -l
+echo height_ngp2_robust_int1_fv_2df.out.txt && diff height_ngp2_robust_int1_fv_2df.out.txt  test_other/height_ngp2_robust_int1_fv_2df.out.txt |wc -l
\ No newline at end of file

checks/inputfiles/README

+This README file describes the 'special' properties of the SNPs in the
+test input files, as well as any other relevant information related
+the files in the inputfiles directory.
+
+* The SNPs
+  - SNP 1: Tests multi-letter alleles (indel, e.g. from 1kG) in Al1 (bug
+           #2586).
+  - SNP 2: Tests multi-letter alleles (indel, e.g. from 1kG) in Al2 (bug
+           #2586); Tests NaN in genetic data.
+  - SNP 3:
+  - SNP 4:
+  - SNP 5: Tests low Rsq in info file.
+  - SNP 6: Tests low MAF, only one out of 200 people has a different dosage.
+
+* Other files

checks/inputfiles/mmscore_InvSigma_aj.sex.age.small.dat

+id4 0.0182686989332702 -1.85665437227337e-05 3.53271468398681e-05 -1.74067659984958e-05 7.43017952446252e-05 -4.50104634347183e-06 -0.000124985264981788 3.3678291247208e-05 5.91968582288439e-05 -0.000243645201364358
+id10 -1.85665437227311e-05 0.0182033225536506 -1.97096076622397e-05 -5.96162563595784e-05 -5.4476919601606e-05 -9.50226926118848e-05 8.62259783546112e-05 7.3858364200391e-05 3.10831766033596e-05 -0.000130544230147021
+id25 3.53271468398634e-05 -1.97096076622424e-05 0.0181688614380818 8.52863460066612e-05 -2.43025930791677e-05 2.81566776967037e-05 4.85600196715389e-06 -6.56752842697857e-06 -1.78033754378254e-05 -1.72662880972169e-05
+id33 -1.74067659984974e-05 -5.96162563595756e-05 8.52863460066646e-05 0.0181194433068409 -5.18411780980829e-05 1.50204770703661e-05 8.56024933212197e-06 4.86573389645691e-05 3.91817984047859e-05 2.79910558167308e-05
+id35 7.43017952446215e-05 -5.44769196016117e-05 -2.43025930791681e-05 -5.18411780980874e-05 0.018106435724767 2.19056691758357e-05 -0.000128600049471664 7.17399863803957e-06 -4.24308810003088e-05 -9.37329781039098e-06
+id58 -4.50104634347325e-06 -9.50226926118893e-05 2.81566776967108e-05 1.50204770703701e-05 2.19056691758332e-05 0.0179881874374435 -1.17277182399346e-05 -4.99197036924978e-05 5.71784745599148e-05 7.8067768163972e-05
+id68 -0.000124985264981788 8.62259783546178e-05 4.85600196715395e-06 8.56024933212903e-06 -0.000128600049471660 -1.17277182399364e-05 0.0181586089581595 2.04445273957100e-05 9.39041676761655e-05 2.69629411670812e-05
+id74 3.36782912472085e-05 7.38583642003905e-05 -6.56752842698657e-06 4.86573389645756e-05 7.17399863804557e-06 -4.99197036925027e-05 2.04445273957019e-05 0.0179962093643653 -0.000143202021910123 -0.000129318612519593
+id107 5.91968582288441e-05 3.10831766033532e-05 -1.78033754378163e-05 3.91817984047789e-05 -4.2430881000305e-05 5.71784745599093e-05 9.3904167676175e-05 -0.000143202021910134 0.0182122458019134 -8.89499732951819e-05
+id119 -0.000243645201364356 -0.000130544230147016 -1.72662880972305e-05 2.79910558167332e-05 -9.37329781039559e-06 7.80677681639793e-05 2.69629411670782e-05 -0.000129318612519587 -8.89499732951819e-05 0.0183196994000799