Commit 1ff22c5e authored by Andreas Tille's avatar Andreas Tille

Update upstream source from tag 'upstream/0.5.0'

Update to upstream version '0.5.0'
with Debian dir 70ecc86c1d02939b4580ba7975fac2f4b3dd9e1e
parents 23a3ee06 495d053f
# /
/.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
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
This diff is collapsed.
* 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
This diff is collapsed.
......@@ -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)
This diff is collapsed.
......@@ -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 \
......
This diff is collapsed.
......@@ -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",
......
......@@ -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
......
......@@ -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")