[Genabel-commits] r1995 - in pkg/ProbABEL: . checks checks/R-tests checks/inputfiles doc src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Jun 16 15:55:40 CEST 2015
Author: lckarssen
Date: 2015-06-16 15:55:39 +0200 (Tue, 16 Jun 2015)
New Revision: 1995
Added:
pkg/ProbABEL/checks/inputfiles/coxph_data_nocovar.txt
pkg/ProbABEL/checks/test_cox_nocovar.sh
Removed:
pkg/ProbABEL/checks/check_pacoxph_nocovar.sh
Modified:
pkg/ProbABEL/checks/Makefile.am
pkg/ProbABEL/checks/R-tests/Makefile.am
pkg/ProbABEL/checks/R-tests/run_model_coxph.R
pkg/ProbABEL/checks/R-tests/run_models_in_R_pacox.R
pkg/ProbABEL/configure.ac
pkg/ProbABEL/doc/ChangeLog
pkg/ProbABEL/doc/ProbABEL_manual.tex
pkg/ProbABEL/doc/pacoxph.1
pkg/ProbABEL/doc/palinear.1
pkg/ProbABEL/doc/palogist.1
pkg/ProbABEL/doc/probabel.1
pkg/ProbABEL/src/command_line_settings.cpp
pkg/ProbABEL/src/coxph_data.cpp
pkg/ProbABEL/src/phedata.cpp
Log:
Merged the ProbABEL v0.4.4-hotfix branch (i.e. ProbABEL v0.4.5) into trunk.
(using: svn merge -r1925:1993 ^/branches//ProbABEL-v0.4.4-hotfix/ProbABEL)
Modified: pkg/ProbABEL/checks/Makefile.am
===================================================================
--- pkg/ProbABEL/checks/Makefile.am 2015-06-16 13:14:53 UTC (rev 1994)
+++ pkg/ProbABEL/checks/Makefile.am 2015-06-16 13:55:39 UTC (rev 1995)
@@ -17,6 +17,7 @@
phenofiles = $(input_files_dir)/height.txt \
$(input_files_dir)/logist_data.txt \
$(input_files_dir)/coxph_data.txt \
+ $(input_files_dir)/coxph_data_nocovar.txt \
$(input_files_dir)/mmscore_pheno.PHE
genofiles = $(input_files_dir)/test.mldose \
@@ -49,7 +50,7 @@
check_SCRIPTS += test_bt.sh
endif
if BUILD_pacoxph
-check_SCRIPTS += test_cox.sh check_pacoxph_nocovar.sh
+check_SCRIPTS += test_cox.sh test_cox_nocovar.sh
endif
AM_TESTS_ENVIRONMENT= \
@@ -133,6 +134,18 @@
coxph_prob_recess.out.txt coxph_prob_add.out.txt \
coxph_prob_2df.out.txt coxph_prob_fv_2df.out.txt
+cleanfiles_cox_nocovar = coxph_dose_nocovar_add.out.txt \
+coxph_dose_nocovar_fv_add.out.txt \
+coxph_prob_nocovar_fv_over_domin.out.txt \
+coxph_prob_nocovar_fv_domin.out.txt \
+coxph_prob_nocovar_over_domin.out.txt \
+coxph_prob_nocovar_fv_add.out.txt \
+coxph_prob_nocovar_fv_recess.out.txt \
+coxph_prob_nocovar_domin.out.txt \
+coxph_prob_nocovar_recess.out.txt \
+coxph_prob_nocovar_add.out.txt \
+coxph_prob_nocovar_2df.out.txt \
+coxph_prob_nocovar_fv_2df.out.txt
dose_files = chr1.chunk1.dose chr1.chunk2.dose chr2.chunk1.dose \
chr2.chunk2.dose chr1.dose chr2.dose
@@ -155,4 +168,5 @@
cleanfiles_dose_check = test.mldose
CLEANFILES = $(cleanfiles_probabel_check) $(cleanfiles_dose_check) \
- $(cleanfiles_bt) $(cleanfiles_qt) $(cleanfiles_mms) $(cleanfiles_cox)
+ $(cleanfiles_bt) $(cleanfiles_qt) $(cleanfiles_mms) $(cleanfiles_cox) \
+ $(cleanfiles_cox_nocovar)
Modified: pkg/ProbABEL/checks/R-tests/Makefile.am
===================================================================
--- pkg/ProbABEL/checks/R-tests/Makefile.am 2015-06-16 13:14:53 UTC (rev 1994)
+++ pkg/ProbABEL/checks/R-tests/Makefile.am 2015-06-16 13:55:39 UTC (rev 1995)
@@ -99,5 +99,18 @@
coxph_prob_recess.out.txt coxph_prob_add.out.txt \
coxph_prob_2df.out.txt coxph_prob_fv_2df.out.txt
+cleanfiles_cox_nocovar = coxph_dose_nocovar_add.out.txt \
+coxph_dose_nocovar_fv_add.out.txt \
+coxph_prob_nocovar_fv_over_domin.out.txt \
+coxph_prob_nocovar_fv_domin.out.txt \
+coxph_prob_nocovar_over_domin.out.txt \
+coxph_prob_nocovar_fv_add.out.txt \
+coxph_prob_nocovar_fv_recess.out.txt \
+coxph_prob_nocovar_domin.out.txt \
+coxph_prob_nocovar_recess.out.txt \
+coxph_prob_nocovar_add.out.txt \
+coxph_prob_nocovar_2df.out.txt \
+coxph_prob_nocovar_fv_2df.out.txt
+
CLEANFILES = $(cleanfiles_bt) $(cleanfiles_qt) $(cleanfiles_mms) \
- $(cleanfiles_cox)
+ $(cleanfiles_cox) $(cleanfiles_cox_nocovar)
Modified: pkg/ProbABEL/checks/R-tests/run_model_coxph.R
===================================================================
--- pkg/ProbABEL/checks/R-tests/run_model_coxph.R 2015-06-16 13:14:53 UTC (rev 1994)
+++ pkg/ProbABEL/checks/R-tests/run_model_coxph.R 2015-06-16 13:55:39 UTC (rev 1995)
@@ -15,8 +15,8 @@
##' output can be compared to the ProbABEL output.
##' @author L.C. Larsen
run.model <- function(model0.txt, model.txt,
- snpcomponent1, snpcomponent2="1") {
-
+ snpcomponent1, snpcomponent2="1",
+ verbose=FALSE) {
if (snpcomponent2 != "1") {
## SNP component 2 is not constant: assume we run the 2df
## model.
@@ -28,15 +28,31 @@
resultR <- data.frame()
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))
snp2 <- eval(parse(text=snpcomponent2))
snp <- snp1 + snp2
+ mu <- rep(1., length(snp)) # Add a constant mean
noNA <- which( !is.na(snp) )
model.0 <- eval(parse(text=model0.txt))
+ ## Check the imputation R^2, if below threshold ProbABEL will
+ ## set the coefficients to NaN and skip the rest of this
+ ## iteration of the for loop.
+ rsq <- Rsq[i-2]
+ if (rsq < rsq.thresh ) {
+ if (twoDF) {
+ row <- c(rsq, NaN, NaN, NaN, NaN, NaN)
+ } else {
+ row <- c(rsq, NaN, NaN, NaN)
+ }
+ resultR <- rbind(resultR, row)
+ next
+ }
+
## Evaluate the model. The whole tryCatch is needed to catch
## problems with non-converging regression.
model = tryCatch({
@@ -49,11 +65,25 @@
eval(parse(text=model.txt)),
war)
)
- })
+ },
+ error = function(err) {
+ return(list(
+ "Can't calculate model; some error occurred",
+ err)
+ )
+ })
- if ( grepl("infinite", model[[2]]$message) |
- grepl("singular", model[[2]]$message) ) {
- ## The model did not converge, fill the coefficients with
+ if (verbose) print(model)
+
+ 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)
+ )
+ ) {
+ ## The model did not converge or some other
+ ## errors/warnings occurred, fill the coefficients with
## NaNs
if (twoDF) {
smA1A2 <- c(NaN, NaN)
@@ -66,30 +96,21 @@
## No convergence problems, we can trust the
## coefficients.
coeff <- summary(model[[1]])$coefficients
+ # The SNP coefficient(s) is/are in the last row(s) of coeff
+ nr <- nrow(coeff)
if (twoDF) {
- smA1A2 <- coeff[4, c("coef", "se(coef)")]
- smA1A1 <- coeff[5, c("coef", "se(coef)")]
+ smA1A2 <- coeff[nr - 1, c("coef", "se(coef)")]
+ smA1A1 <- coeff[nr, c("coef", "se(coef)")]
} else {
- sm <- coeff[4, c("coef", "se(coef)")]
+ sm <- coeff[nr, c("coef", "se(coef)")]
}
lrt <- 2 * ( model[[1]]$loglik[2] - model.0$loglik[2] )
}
- ## Check the imputation R^2, if below threshold ProbABEL will
- ## set the coefficients to NaN.
- rsq <- Rsq[i-2]
if (twoDF) {
- if( rsq < rsq.thresh ) {
- row <- c(rsq, NaN, NaN, NaN, NaN, NaN)
- } else {
- row <- c(rsq, smA1A2[1], smA1A2[2], smA1A1[1], smA1A1[2], lrt)
- }
+ row <- c(rsq, smA1A2[1], smA1A2[2], smA1A1[1], smA1A1[2], lrt)
} else {
- if( rsq < rsq.thresh ) {
- row <- c(rsq, NaN, NaN, NaN)
- } else {
- row <- c(rsq, sm[1], sm[2], lrt)
- }
+ row <- c(rsq, sm[1], sm[2], lrt)
}
resultR <- rbind(resultR, row)
Modified: pkg/ProbABEL/checks/R-tests/run_models_in_R_pacox.R
===================================================================
--- pkg/ProbABEL/checks/R-tests/run_models_in_R_pacox.R 2015-06-16 13:14:53 UTC (rev 1994)
+++ pkg/ProbABEL/checks/R-tests/run_models_in_R_pacox.R 2015-06-16 13:55:39 UTC (rev 1995)
@@ -53,8 +53,8 @@
source(paste0(srcdir, "run_model_coxph.R"))
model.fn.0 <-
- "coxph( Surv(fupt_chd, chd)[noNA] ~ sex[noNA] + age[noNA] + othercov[noNA] )"
-model.fn <- "coxph( Surv(fupt_chd, chd) ~ sex + age + othercov + snp1 )"
+ "coxph( Surv(fupt_chd, chd)[noNA] ~ mu[noNA] + sex[noNA] + age[noNA] + othercov[noNA] )"
+model.fn <- "coxph( Surv(fupt_chd, chd) ~ mu + sex + age + othercov + snp1 )"
## Additive model, dosages
prnt(" additive (dosages)")
@@ -112,3 +112,90 @@
rownames(prob.2df.R) <- NULL
stopifnot( all.equal(prob.2df.PA, prob.2df.R, tol=tol) )
cat("OK\n")
+
+
+
+cat("\nRun the same checks, but without any covariates (see bug #1266)\n")
+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)[, colsAdd]
+prob.add.PA <- read.table("coxph_prob_nocovar_add.out.txt",
+ 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",
+ head=TRUE)[, colsRec]
+prob.odom.PA <- read.table("coxph_prob_nocovar_over_domin.out.txt",
+ head=TRUE)[, colsOdom]
+prob.2df.PA <- read.table("coxph_prob_nocovar_2df.out.txt",
+ head=TRUE)[, cols2df]
+
+model.fn.0 <-
+ "coxph( Surv(fupt_chd, chd)[noNA] ~ mu[noNA] )"
+model.fn <- "coxph( Surv(fupt_chd, chd) ~ mu + snp1 )"
+
+
+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) <- colsAdd
+rownames(dose.add.R) <- NULL
+stopifnot( all.equal(dose.add.PA, dose.add.R, tol=tol) )
+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) <- colsAdd
+rownames(prob.add.R) <- NULL
+stopifnot( all.equal(prob.add.PA, prob.add.R, tol=tol) )
+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("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("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("OK\n")
+
+
+## 2df model
+prnt(" 2df")
+model.fn <-
+ "coxph( Surv(fupt_chd, chd) ~ mu+ snp1 + snp2 )"
+snpd1 <- "prob[, indexHet]"
+snpd2 <- "prob[, indexHom]"
+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("OK\n")
Deleted: pkg/ProbABEL/checks/check_pacoxph_nocovar.sh
===================================================================
--- pkg/ProbABEL/checks/check_pacoxph_nocovar.sh 2015-06-16 13:14:53 UTC (rev 1994)
+++ pkg/ProbABEL/checks/check_pacoxph_nocovar.sh 2015-06-16 13:55:39 UTC (rev 1995)
@@ -1,50 +0,0 @@
-#!/bin/bash
-# L.C. Karssen
-# This script is used to test whether the Cox PH regression works
-# correctly when no covariate is present in the phenotype input file.
-# Currently this test fails, see bug #1266.
-
-echo "Testing pacoxph without covariates..."
-
-if [ -z ${AWK} ]; then
- AWK=awk
-fi
-
-# Exit with error when one of the steps in the script fails
-set -e
-
-# -------- Set some default paths and file names -------
-if [ -z ${srcdir} ]; then
- srcdir="."
-fi
-inputdir="${srcdir}/inputfiles/"
-padir="${srcdir}/../src/"
-
-dosefile="$inputdir/test.mldose"
-infofile="$inputdir/test.mlinfo"
-mapfile="$inputdir/test.map"
-orig_phenofile="$inputdir/coxph_data.txt"
-phenofile="coxph_data.txt"
-outfile="pacoxph_nocovar"
-pacoxph="${padir}/pacoxph"
-
-# ------ Prepare the phenotype file by removing the covariate column
-# from the existing phenotype file ------
-${AWK} '{print $1, $2, $3}' $orig_phenofile > $phenofile
-
-
-# ---------- function definitions ----------
-run_test ()
-{
- ## When bug #1266 is fixed, this function should be expanded to
- ## include a verification against known-good results.
- echo "Checking whether Cox PH regression works without covariates..."
- $pacoxph -p $phenofile -d $dosefile -i $infofile -m $mapfile \
- -o $outfile
-
-}
-
-
-# ---------- Continuation of the main function ----------
-run_test
-echo "---- Finished check for Cox regression without covariates ----"
Copied: pkg/ProbABEL/checks/inputfiles/coxph_data_nocovar.txt (from rev 1993, branches/ProbABEL-v0.4.4-hotfix/ProbABEL/checks/inputfiles/coxph_data_nocovar.txt)
===================================================================
--- pkg/ProbABEL/checks/inputfiles/coxph_data_nocovar.txt (rev 0)
+++ pkg/ProbABEL/checks/inputfiles/coxph_data_nocovar.txt 2015-06-16 13:55:39 UTC (rev 1995)
@@ -0,0 +1,201 @@
+id fupt_chd chd
+id636728 3.187930645 0
+id890314 2.099691952 0
+id102874 9.133488079 1
+id200949 7.525406804 0
+id336491 6.798229522 0
+id988766 6.149545358 0
+id21999 1.013546103 1
+id433893 1.282853098 0
+id688932 8.340206657 0
+id394203 3.392345681 1
+id995678 9.587748608 0
+id694339 1.25827049 0
+id256455 5.414840139 0
+id14836 6.177914833 1
+id817128 2.940033056 1
+id803325 9.619267566 0
+id521287 0.320487323 1
+id701472 7.101945547 0
+id850010 6.780876558 1
+id268483 7.903947419 0
+id738781 1.086018804 0
+id28411 1.254350142 1
+id541635 9.057035035 1
+id751101 1.303018406 0
+id826300 6.735484141 1
+id884387 9.097303962 0
+id492414 4.756845946 0
+id268871 4.018189958 0
+id627354 8.874035732 0
+id503932 0.091310535 0
+id163442 3.240808272 0
+id317797 9.050109718 0
+id687857 7.992764162 0
+id871570 2.353719529 0
+id724067 1.617985994 1
+id874076 3.099430273 0
+id927863 3.64738684 0
+id369805 0.973580503 0
+id668376 0.567035975 0
+id717362 6.77104712 1
+id665504 8.956792361 0
+id336637 7.196501733 0
+id60633 5.086557706 0
+id848600 8.270336349 1
+id169514 9.421256192 0
+id690732 9.952887485 0
+id684760 3.361963101 0
+id553502 5.557606017 1
+id214917 6.535326783 0
+id849169 7.903922417 1
+id941921 2.305853849 0
+id784646 3.120883065 1
+id520954 8.756848347 1
+id996355 7.326814334 0
+id96730 6.664821662 0
+id673442 3.081616083 0
+id68305 3.525197877 1
+id653025 7.231586593 0
+id208543 4.622250109 0
+id335725 5.924295793 1
+id980400 2.35961827 0
+id869939 0.736202662 0
+id297563 6.641749611 0
+id852663 5.124175591 0
+id162070 7.910906314 0
+id272875 9.825455416 1
+id163787 2.471002157 0
+id422204 2.915218179 0
+id120197 3.367355865 0
+id33660 4.407010904 1
+id803855 8.755961077 0
+id255048 7.749378138 0
+id690936 7.543593942 0
+id126807 0.251036587 0
+id99016 6.491002725 0
+id883847 1.088929613 0
+id354523 9.618060855 1
+id737255 4.133744506 1
+id990941 5.646270286 0
+id25464 8.272429793 0
+id918375 0.723567779 0
+id537828 4.0147246 0
+id682778 9.130628394 0
+id587547 7.44191196 0
+id670874 4.512303089 0
+id444459 8.723388559 0
+id777456 4.115325335 1
+id452384 2.067203142 0
+id826975 8.348432332 0
+id519567 6.785076383 0
+id84292 7.634117891 0
+id124432 5.099839207 1
+id800145 8.551368147 0
+id153857 3.556301609 0
+id587157 8.758929135 1
+id506262 6.238330685 1
+id634462 0.905623161 0
+id687592 0.290071515 0
+id955526 1.838730964 0
+id181850 0.144560173 1
+id159506 6.045008786 1
+id609051 6.873427527 1
+id963886 0.853758289 0
+id405792 3.909860553 1
+id494172 7.480072965 0
+id964637 3.816858446 0
+id799355 1.226117869 1
+id157111 4.166491482 1
+id114524 8.645968152 0
+id954931 0.158490245 0
+id827034 0.702377072 0
+id689645 2.0658374 0
+id281585 4.640052191 1
+id885624 3.615978697 0
+id577871 2.230068497 0
+id238796 5.161400753 0
+id481035 7.693009686 0
+id972713 6.34411768 1
+id905484 4.55801149 0
+id713511 6.721560405 1
+id512328 4.431789868 0
+id703534 8.974882244 1
+id409904 4.845419164 0
+id577169 1.384047668 0
+id813971 6.785385478 0
+id558483 0.27946972 0
+id892784 1.505769901 0
+id611178 3.97305697 1
+id192732 7.859877791 0
+id917280 1.39759155 0
+id435876 4.928142177 0
+id980722 0.72730016 1
+id308273 3.466919258 1
+id476685 8.925909139 0
+id315883 7.79623385 0
+id935945 3.491058239 1
+id991781 6.268341084 0
+id65199 2.65595714 0
+id226233 6.409789501 0
+id860183 7.567894972 0
+id295209 0.983143735 1
+id544964 1.927387219 0
+id648663 7.914653713 0
+id710165 9.981244413 0
+id392593 7.418579368 0
+id129945 9.344802496 0
+id382621 1.828452452 0
+id901440 9.496967232 0
+id39847 8.962161053 0
+id526460 1.855905599 0
+id477473 1.89081994 1
+id448194 4.814151059 0
+id904184 8.028988991 1
+id747852 3.836260637 1
+id711012 0.341281233 0
+id683879 9.297243459 0
+id789575 5.098990324 0
+id650729 7.941450537 0
+id934302 9.134082902 1
+id555013 8.812444844 1
+id82779 7.657743527 0
+id771444 6.664855358 1
+id821562 5.19036265 0
+id292809 9.909691181 1
+id645690 6.515405399 0
+id223901 6.355198262 1
+id41320 0.397540162 1
+id96181 3.155308448 1
+id147900 4.719584771 0
+id702917 7.72554794 0
+id150640 6.511110099 1
+id518391 9.626758766 0
+id879076 7.61506094 0
+id952031 4.87888217 1
+id10055 5.724749579 0
+id727213 5.731963291 1
+id41961 6.818411434 0
+id257209 9.897638239 1
+id995361 3.151184642 0
+id957918 5.30099154 0
+id975370 9.047264812 0
+id889896 9.578142144 1
+id978164 4.726348119 0
+id90359 6.237920077 0
+id307158 5.49680528 0
+id755940 1.683472511 0
+id995582 1.486566562 0
+id363965 7.723182293 1
+id729124 6.861379688 0
+id871963 1.399526332 1
+id475172 2.745065294 0
+id804699 8.378029915 1
+id625843 4.925798265 1
+id595713 3.273883772 0
+id462604 5.541305429 0
+id106141 9.529928524 0
+id689349 9.243194854 0
+id639003 6.594759182 1
+id393896 8.69916317 1
+id450307 7.291170542 1
Copied: pkg/ProbABEL/checks/test_cox_nocovar.sh (from rev 1993, branches/ProbABEL-v0.4.4-hotfix/ProbABEL/checks/test_cox_nocovar.sh)
===================================================================
--- pkg/ProbABEL/checks/test_cox_nocovar.sh (rev 0)
+++ pkg/ProbABEL/checks/test_cox_nocovar.sh 2015-06-16 13:55:39 UTC (rev 1995)
@@ -0,0 +1,82 @@
+#!/bin/bash
+# This script runs checks on ProbABEL's pacoxph module
+# These are almost identical to the checks in test_cox.sh, but this
+# time without covariate data (see bug #1266).
+
+
+echo "Analysing Cox model without covariates..."
+
+scriptdir=$(dirname $0)
+
+if [ -z ${PA_BINDIR} ]; then
+ PA_BINDIR="${scriptdir}/../src/"
+fi
+if [ -z ${srcdir} ]; then
+ srcdir="."
+ PA_BINDIR=${scriptdir}/../src/
+fi
+
+. ${scriptdir}/run_diff.sh
+
+inputdir=${scriptdir}/inputfiles
+pacoxph=${PA_BINDIR}/pacoxph
+
+# Redirect all output to file descriptor 3 to /dev/null except if
+# the first argument is "verbose" then redirect handle 3 to stdout
+exec 3>/dev/null
+if [ "$1" = "verbose" ]; then
+ echo "Verbose mode ON"
+ exec 3>&1
+fi
+
+$pacoxph \
+ -p ${inputdir}/coxph_data_nocovar.txt \
+ -d ${inputdir}/test.mldose \
+ -i ${inputdir}/test.mlinfo \
+ -m ${inputdir}/test.map \
+ -c 19 \
+ -o coxph_dose_nocovar \
+ >& 3
+
+$pacoxph \
+ -p ${inputdir}/coxph_data_nocovar.txt \
+ -d ${inputdir}/test.dose.fvi \
+ -i ${inputdir}/test.mlinfo \
+ -m ${inputdir}/test.map \
+ -c 19 \
+ -o coxph_dose_nocovar_fv \
+ >& 3
+
+run_diff coxph_dose_nocovar_add.out.txt coxph_dose_nocovar_fv_add.out.txt \
+ "pacoxph check: dose vs. dose_fv"
+
+
+$pacoxph \
+ -p ${inputdir}/coxph_data_nocovar.txt \
+ -d ${inputdir}/test.mlprob \
+ -i ${inputdir}/test.mlinfo \
+ -m ${inputdir}/test.map \
+ --ngpreds=2 \
+ -c 19 \
+ -o coxph_prob_nocovar \
+ >& 3
+
+run_diff coxph_dose_nocovar_add.out.txt coxph_prob_nocovar_add.out.txt \
+ "pacoxph check: dose vs. prob" -I SNP
+
+
+$pacoxph \
+ -p ${inputdir}/coxph_data_nocovar.txt \
+ -d ${inputdir}/test.prob.fvi \
+ -i ${inputdir}/test.mlinfo \
+ -m ${inputdir}/test.map \
+ --ngpreds=2 \
+ -c 19 \
+ -o coxph_prob_nocovar_fv \
+ >& 3
+
+for model in add domin recess over_domin 2df; do
+ run_diff coxph_prob_nocovar_${model}.out.txt \
+ coxph_prob_nocovar_fv_${model}.out.txt \
+ "pacoxph check ($model model): prob vs. prob_fv"
+done
Modified: pkg/ProbABEL/configure.ac
===================================================================
--- pkg/ProbABEL/configure.ac 2015-06-16 13:14:53 UTC (rev 1994)
+++ pkg/ProbABEL/configure.ac 2015-06-16 13:55:39 UTC (rev 1995)
@@ -1,7 +1,7 @@
# Process this file with autoconf to produce a configure script.
AC_PREREQ([2.67])
-AC_INIT(ProbABEL, 0.4.4, genabel-devel at r-forge.wu-wien.ac.at)
+AC_INIT(ProbABEL, 0.4.5, genabel-devel at r-forge.wu-wien.ac.at)
AM_INIT_AUTOMAKE([silent-rules subdir-objects])
AM_SILENT_RULES([yes])
AC_CONFIG_SRCDIR([src/mlinfo.h])
Modified: pkg/ProbABEL/doc/ChangeLog
===================================================================
--- pkg/ProbABEL/doc/ChangeLog 2015-06-16 13:14:53 UTC (rev 1994)
+++ pkg/ProbABEL/doc/ChangeLog 2015-06-16 13:55:39 UTC (rev 1995)
@@ -1,4 +1,4 @@
-***** v.0.5.0 (2014.)
+***** v.0.5.0 (2015.)
* ProbABEL now depends on the Eigen library. The option to compile without
Eigen has been removed.
* Not-a-number values in the output are now printed as "NaN" instead of
@@ -19,6 +19,26 @@
with sed error. Thanks to forum user mmold for reporting the bug.
+***** v.0.4.5 (2015.05.26)
+* Fixed bug #6054: "Not all ProbABEL's short command line options are
+ correctly parsed."
+* Fixed bug #6041: "ProbABEL's Cox module reports too many errors (beta
+ may be infinite, setting beta and se to 'NaN')". Thanks to Anne
+ Grotenhuis from the Radboud Medical Centre Nijmegen and forum user
+ quentin300 for their extensive bug reports and help in testing the fix.
+* Fixed bug #1266: "pacoxph with no covariates"; ProbABEL now also works
+ when no covariates have been provided. Thanks to Aaron Joon for
+ reporting this bug back in 2011 and to Anne Grotenhuis and to forum user
+ quentin300 for pushing this bug to my attention and their help in beta
+ testing.
+* pacoxph now displays the regression equation correctly.
+* Minor fix to the manual: the Rsq values in the info file should be >
+ 1e-16 (and not > 0 as mentioned before) otherwise the output will be set
+ to 'nan'.
+* Minor fix to the man-pages: according to the man-pages palogist,
+ palinear and pacoxph all did regression using a linear model.
+
+
***** v.0.4.4 (2014.11.07)
* Fixed bug #5729 in the Cox PH module. Some checks for problems with the
regression were incorrectly implemented. Thanks to Matthias Wuttke from
Modified: pkg/ProbABEL/doc/ProbABEL_manual.tex
===================================================================
--- pkg/ProbABEL/doc/ProbABEL_manual.tex 2015-06-16 13:14:53 UTC (rev 1994)
+++ pkg/ProbABEL/doc/ProbABEL_manual.tex 2015-06-16 13:55:39 UTC (rev 1995)
@@ -1,17 +1,17 @@
\documentclass[12pt,a4paper]{article}
-\title{Manual for ProbABEL v0.4.4}
+\title{Manual for ProbABEL v0.4.5}
\author{\emph{Current Programmers:} Lennart Karssen$^{1,2}$, Maarten
Kooyman$^2$, \\
Yurii Aulchenko$^{1,3}$ \\
\emph{Former Programmers:} Maksim Struchalin
\\
\\
- $^{1}${\small PolyOmica} \\
- $^{2}${\small Erasmus MC, Rotterdam}\\
+ $^{1}${\small PolyOmica, Groningen, The Netherlands} \\
+ $^{2}${\small Erasmus MC, Rotterdam, The Netherlands}\\
$^{3}${\small Institute of Cytology and Genetics SD RAS, Novosibirsk}
}
-\date{November 7, 2014}
+\date{May 26, 2015}
\usepackage[utf8]{inputenc}
@@ -221,15 +221,15 @@
multiplication. In ProbABEL versions before v0.5.0 EIGEN was not
required (but still recommended).}. In order to download the
library, go to \url{http://eigen.tuxfamily.org} and download the
-\texttt{tar.gz} file of the latest version of EIGEN (3.2.1 at the time
+\texttt{tar.gz} file of the latest version of EIGEN (3.2.4 at the time
of writing). Extract the files:
\begin{lstlisting}[basicstyle=\footnotesize\ttfamily,]
-user at server:~$ tar -xzf 3.2.1.tar.gz
+user at server:~$ tar -xzf 3.2.4.tar.gz
\end{lstlisting}
This will create a directory called \texttt{eigen-eigen} followed by a
series of letters and digits. For simplicity we rename it to EIGEN
\begin{lstlisting}[basicstyle=\footnotesize\ttfamily,]
-user at server:~$ mv eigen-eigen-6b38706d90a9 EIGEN
+user at server:~$ mv eigen-eigen-10219c95fe65 EIGEN
\end{lstlisting}
Now it's time to extract the \PA{} source code and move into the
@@ -304,12 +304,13 @@
Actually, for \PA{}, it (almost) does not matter what is written in
this file -- this information is simply copied to the output. However,
\textbf{it is critical} that the number of columns is
-seven\footnote{This means that for \texttt{minimac} output files the number of
- columns needs to be reduced. This can be done using e.g.~GAWK or
- \texttt{cut}.} and the number of lines in the file is equal to the
-number of SNPs in the corresponding DOSE file (plus one for the header
-line). Also make sure that the ``Rsq'' column contains values $>0$,
-otherwise you will end up with $\beta$'s set to \texttt{nan}.
+seven\footnote{This means that for \texttt{minimac} output files the
+ number of columns needs to be reduced. This can be done using
+ e.g.~GAWK or \texttt{cut}.} and the number of lines in the file is
+equal to the number of SNPs in the corresponding DOSE file (plus one
+for the header line). Also make sure that the ``Rsq'' column contains
+values $>1 \cdot 10^{-16}$, otherwise you will end up with $\beta$'s
+set to \texttt{nan}.
The example of SNP information file content follows here (also
to be found in \texttt{ProbABEL/examples/test.mlinfo})
@@ -461,7 +462,7 @@
short explanation to the command line options:
\begin{verbatim}
user at server:~$ palogist --help
-probabel v. 0.4.4
+probabel v. 0.4.5
(C) Yurii Aulchenko, Lennart C. Karssen, Maksim Struchalin, EMCR
Using EIGEN version 3.1.2 for matrix operations
@@ -884,20 +885,20 @@
Abecasis (2007). The general analysis model is a linear mixed model
where the expectation of the trait is defined as
$$
-E[\mathbf{Y}] = \mathbf{X} \mathbf{\beta},
+E[\mathbf{Y}] = \mathbf{X} \boldsymbol{\beta},
$$
identical to that defined for a linear model
(cf.~Eq.~\ref{eq:expectation}). To account for correlations between
the phenotypes of relatives which may be induced by family relations
the variance-covariance matrix is defined to be proportional to the
linear combination of the identity matrix $\mathbf{I}$ and the
-relationship matrix $\mathbf{\Phi}$:
+relationship matrix $\boldsymbol{\Phi}$:
$$
-\mathbf{V}_{\sigma^2,h^2} = \sigma^2 \left( 2 h^2 \mathbf{\Phi} + (1-h^2)
+\mathbf{V}_{\sigma^2,h^2} = \sigma^2 \left( 2 h^2 \boldsymbol{\Phi} + (1-h^2)
\mathbf{I} \right),
$$
where $h^2$ is the heritability of the trait. The relationship matrix
-$\mathbf{\Phi}$ is twice the matrix containing the coefficients of
+$\boldsymbol{\Phi}$ is twice the matrix containing the coefficients of
kinship between all pairs of individuals under consideration; its
estimation is discussed separately in section \ref{kinship}.
@@ -916,19 +917,19 @@
models fit in GWAS (e.g.\ effects of sex, age, etc.), and the part
including SNP information, $\mathbf{X_g}$:
$$
-E[\mathbf{Y}] = \mathbf{X}_x \mathbf{\beta}_x +
-\mathbf{X}_g \mathbf{\beta}_g.
+E[\mathbf{Y}] = \mathbf{X}_x \boldsymbol{\beta}_x +
+\mathbf{X}_g \boldsymbol{\beta}_g.
$$
Note that the latter design matrix may include not only the main SNP
effect, but e.g.\ SNP by environment interaction terms.
In the first step, a linear mixed model not including SNP effects
$$
-E[\mathbf{Y}] = \mathbf{X}_x \mathbf{\beta}_x
+E[\mathbf{Y}] = \mathbf{X}_x \boldsymbol{\beta}_x
$$
is fitted. The maximum likelihood estimates (MLEs) of the model
parameters (regression coefficients for the fixed effects
-$\hat{\mathbf{\beta}}_x$, the residual variance $\hat{\sigma}^2_x$ and
+$\hat{\boldsymbol{\beta}}_x$, the residual variance $\hat{\sigma}^2_x$ and
the heritability $\hat{h}^2_x$) can be obtained by numerical
maximization of the likelihood function
$$
@@ -984,7 +985,7 @@
\subsubsection{Estimation of the kinship matrix}
\label{kinship}
-The relationship matrix $\mathbf{\Phi}$ used in estimation of the
+The relationship matrix $\boldsymbol{\Phi}$ used in estimation of the
linear mixed model for pedigree data is twice the matrix containing
the coefficients of kinship between all pairs of individuals under consideration.
This coefficient is defined as the probability that two gametes randomly sampled
@@ -1035,7 +1036,7 @@
\end{quote}
A proper reference may look like
\begin{quote}
-For the analysis of imputed data, we used \PA{} v.0.4.4
+For the analysis of imputed data, we used \PA{} v.0.4.5
from the \texttt{GenABEL} suite of programs (Aulchenko \emph{et al.}, 2010).
\end{quote}
@@ -1080,3 +1081,8 @@
\printindex
\end{document}
+
+%%% Local Variables:
+%%% mode: latex
+%%% TeX-master: t
+%%% End:
Modified: pkg/ProbABEL/doc/pacoxph.1
===================================================================
--- pkg/ProbABEL/doc/pacoxph.1 2015-06-16 13:14:53 UTC (rev 1994)
+++ pkg/ProbABEL/doc/pacoxph.1 2015-06-16 13:55:39 UTC (rev 1995)
@@ -1,6 +1,7 @@
-.TH pacoxph 1 "7 November 2014" "ProbABEL 0.4.4"
+.TH pacoxph 1 "26 May 2015" "ProbABEL 0.4.5"
.SH NAME
-pacoxph \- Perform Genome-Wide Association Analysis using a linear model
+pacoxph \- Perform Genome-Wide Association Analysis using Cox' Proportional
+hazards model
.SH SYNOPSIS
.B pacoxph
.RI "[ " "command-line options" " ]"
Modified: pkg/ProbABEL/doc/palinear.1
===================================================================
--- pkg/ProbABEL/doc/palinear.1 2015-06-16 13:14:53 UTC (rev 1994)
+++ pkg/ProbABEL/doc/palinear.1 2015-06-16 13:55:39 UTC (rev 1995)
@@ -1,4 +1,4 @@
-.TH palinear 1 "7 November 2014" "ProbABEL 0.4.4"
+.TH palinear 1 "26 May 2015" "ProbABEL 0.4.5"
.SH NAME
palinear \- Perform Genome-Wide Association Analysis using a linear model
.SH SYNOPSIS
Modified: pkg/ProbABEL/doc/palogist.1
===================================================================
--- pkg/ProbABEL/doc/palogist.1 2015-06-16 13:14:53 UTC (rev 1994)
+++ pkg/ProbABEL/doc/palogist.1 2015-06-16 13:55:39 UTC (rev 1995)
@@ -1,6 +1,6 @@
-.TH palogist 1 "7 November 2014" "ProbABEL 0.4.4"
+.TH palogist 1 "26 May 2015" "ProbABEL 0.4.5"
.SH NAME
-palogist \- Perform Genome-Wide Association Analysis using a linear model
+palogist \- Perform Genome-Wide Association Analysis using a logistic model
.SH SYNOPSIS
.B palogist
.RI "[ " "command-line options" " ]"
Modified: pkg/ProbABEL/doc/probabel.1
===================================================================
--- pkg/ProbABEL/doc/probabel.1 2015-06-16 13:14:53 UTC (rev 1994)
+++ pkg/ProbABEL/doc/probabel.1 2015-06-16 13:55:39 UTC (rev 1995)
@@ -1,4 +1,4 @@
-.TH ProbABEL 1 "7 November 2014" "ProbABEL 0.4.4"
+.TH ProbABEL 1 "26 May 2015" "ProbABEL 0.4.5"
.SH NAME
probabel \- Wrapper around the three ProbABEL binaries, simplifying their use
.SH SYNOPSIS
Modified: pkg/ProbABEL/src/command_line_settings.cpp
===================================================================
--- pkg/ProbABEL/src/command_line_settings.cpp 2015-06-16 13:14:53 UTC (rev 1994)
+++ pkg/ProbABEL/src/command_line_settings.cpp 2015-06-16 13:55:39 UTC (rev 1995)
@@ -159,7 +159,7 @@
void cmdvars::set_variables(int argc, char * argv[])
{
int next_option;
- const char * const short_options = "p:i:d:m:n:c:o:s:t:g:a:erlhb:vu";
+ const char * const short_options = "p:i:d:m:n:c:o:s:t:g:a:relhb:k:v:u";
// b - interaction parameter
// ADD --fv FLAG (FILEVECTOR), IN WHICH CASE USE ALTERNATIVE
// CONSTRUCTOR FOR GENDATA
Modified: pkg/ProbABEL/src/coxph_data.cpp
===================================================================
--- pkg/ProbABEL/src/coxph_data.cpp 2015-06-16 13:14:53 UTC (rev 1994)
+++ pkg/ProbABEL/src/coxph_data.cpp 2015-06-16 13:55:39 UTC (rev 1995)
@@ -143,13 +143,14 @@
exit(1);
}
- X.reinit(nids, ncov);
+ X.reinit(nids, (ncov + 1));
stime.reinit(nids, 1);
sstat.reinit(nids, 1);
weights.reinit(nids, 1);
offset.reinit(nids, 1);
strata.reinit(nids, 1);
order.reinit(nids, 1);
+
for (int i = 0; i < nids; i++)
{
stime[i] = (phed.Y).get(i, 0);
@@ -163,12 +164,21 @@
}
}
- for (int j = 0; j < phed.ncov; j++)
+ // Add a column with a constant (=1) to the X matrix (the mean)
+ for (int i = 0; i < nids; i++)
{
+ X.put(1., i, 0);
+ }
+
+ for (int j = 1; j <= phed.ncov; j++)
+ {
for (int i = 0; i < nids; i++)
- X.put((phed.X).get(i, j), i, j);
+ {
+ X.put((phed.X).get(i, j - 1), i, j);
+ }
}
+
if (snpnum > 0)
{
for (int j = 0; j < ngpreds; j++)
@@ -285,7 +295,7 @@
gend->get_var(snpnum * ngpreds + j, snpdata);
for (int i = 0; i < nids; i++) {
- X.put(snpdata[i], (ncov - j - 1), order[i]);
+ X.put(snpdata[i], (ncov - j), order[i]);
if (std::isnan(snpdata[i])) {
masked_data[order[i]] = true;
// snp not masked
@@ -475,12 +485,25 @@
// to NA if flag < nvar (with nvar = ncol(x)) and MAXITER >
// 0. coxph.R then checks for any NAs in the coefficients and
// outputs the warning message if NAs were found.
- if (flag < X.nrow && MAXITER > 0) {
- cerr << "Warning for " << snpinfo.name[cursnp]
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/genabel -r 1995
More information about the Genabel-commits
mailing list