[Genabel-commits] r1682 - in pkg/ProbABEL: . checks checks/R-tests checks/verified_results examples src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Apr 7 18:34:40 CEST 2014
Author: lckarssen
Date: 2014-04-07 18:34:40 +0200 (Mon, 07 Apr 2014)
New Revision: 1682
Modified:
pkg/ProbABEL/checks/R-tests/initial_checks.R
pkg/ProbABEL/checks/R-tests/run_models_in_R_pacox.R
pkg/ProbABEL/checks/R-tests/run_models_in_R_palinear.R
pkg/ProbABEL/checks/R-tests/run_models_in_R_palogist.R
pkg/ProbABEL/checks/run_diff.sh
pkg/ProbABEL/checks/test_bt.sh
pkg/ProbABEL/checks/test_cox.sh
pkg/ProbABEL/checks/test_mms.sh
pkg/ProbABEL/checks/test_qt.sh
pkg/ProbABEL/checks/verified_results/height_add.out.txt
pkg/ProbABEL/checks/verified_results/height_base_add.out.txt
pkg/ProbABEL/checks/verified_results/height_ngp2_add.out.txt
pkg/ProbABEL/checks/verified_results/height_ngp2_domin.out.txt
pkg/ProbABEL/checks/verified_results/height_ngp2_over_domin.out.txt
pkg/ProbABEL/checks/verified_results/height_ngp2_recess.out.txt
pkg/ProbABEL/configure.ac
pkg/ProbABEL/examples/mmscore.R
pkg/ProbABEL/src/cholesky.cpp
pkg/ProbABEL/src/command_line_settings.cpp
pkg/ProbABEL/src/command_line_settings.h
pkg/ProbABEL/src/coxph_data.cpp
pkg/ProbABEL/src/coxph_data.h
pkg/ProbABEL/src/data.cpp
pkg/ProbABEL/src/data.h
pkg/ProbABEL/src/fvlib
pkg/ProbABEL/src/gendata.cpp
pkg/ProbABEL/src/gendata.h
pkg/ProbABEL/src/main.cpp
pkg/ProbABEL/src/main_functions_dump.cpp
pkg/ProbABEL/src/maskedmatrix.cpp
pkg/ProbABEL/src/maskedmatrix.h
pkg/ProbABEL/src/probabel
pkg/ProbABEL/src/reg1.cpp
pkg/ProbABEL/src/reg1.h
pkg/ProbABEL/src/regdata.cpp
pkg/ProbABEL/src/regdata.h
Log:
Merged the ProbABEL-0.50 (r1673) branch into trunk.
Modified: pkg/ProbABEL/checks/R-tests/initial_checks.R
===================================================================
--- pkg/ProbABEL/checks/R-tests/initial_checks.R 2014-04-07 15:44:06 UTC (rev 1681)
+++ pkg/ProbABEL/checks/R-tests/initial_checks.R 2014-04-07 16:34:40 UTC (rev 1682)
@@ -65,26 +65,22 @@
## 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",
Modified: pkg/ProbABEL/checks/R-tests/run_models_in_R_pacox.R
===================================================================
--- pkg/ProbABEL/checks/R-tests/run_models_in_R_pacox.R 2014-04-07 15:44:06 UTC (rev 1681)
+++ pkg/ProbABEL/checks/R-tests/run_models_in_R_pacox.R 2014-04-07 16:34:40 UTC (rev 1682)
@@ -31,9 +31,9 @@
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",
@@ -59,7 +59,7 @@
## Additive model, 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 ")
@@ -68,7 +68,7 @@
## Additive model, 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 ")
Modified: pkg/ProbABEL/checks/R-tests/run_models_in_R_palinear.R
===================================================================
--- pkg/ProbABEL/checks/R-tests/run_models_in_R_palinear.R 2014-04-07 15:44:06 UTC (rev 1681)
+++ pkg/ProbABEL/checks/R-tests/run_models_in_R_palinear.R 2014-04-07 16:34:40 UTC (rev 1682)
@@ -20,9 +20,9 @@
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,7 +36,11 @@
## (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)
+#for 2df model the last SNP is interchangeable: EIGEN calculates the beta for the other SNP than R. This causes the beta to have the wrong sign. This part of change the position of the snp beta(and swaps sign) and SE if beta and 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
@@ -53,7 +57,7 @@
## Additive model, 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 ")
@@ -62,7 +66,7 @@
## Additive model, 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 ")
@@ -123,6 +127,7 @@
}
colnames(prob.2df.R) <- cols2df
rownames(prob.2df.R) <- NULL
+
stopifnot( all.equal(prob.2df.PA, prob.2df.R, tol=tol) )
cat("2df\n")
Modified: pkg/ProbABEL/checks/R-tests/run_models_in_R_palogist.R
===================================================================
--- pkg/ProbABEL/checks/R-tests/run_models_in_R_palogist.R 2014-04-07 15:44:06 UTC (rev 1681)
+++ pkg/ProbABEL/checks/R-tests/run_models_in_R_palogist.R 2014-04-07 16:34:40 UTC (rev 1682)
@@ -20,9 +20,9 @@
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",
@@ -52,7 +52,7 @@
## Additive model, 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 ")
@@ -60,7 +60,7 @@
## Additive model, 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 ")
Modified: pkg/ProbABEL/checks/run_diff.sh
===================================================================
--- pkg/ProbABEL/checks/run_diff.sh 2014-04-07 15:44:06 UTC (rev 1681)
+++ pkg/ProbABEL/checks/run_diff.sh 2014-04-07 16:34:40 UTC (rev 1682)
@@ -8,19 +8,13 @@
# $1: first file to compare
# $2: second file to compare
# $3: start message to print on the output line before OK or FAILED
- # $4-: any arguments to pass to the diff command
file1=$1
file2=$2
- name=""
- if [ ${#} -ge 3 ]; then
- name=$3
- shift 3
- args=$@
- fi
+ name=$3
blanks=" "
- if diff $args "$file1" "$file2"; then
+ if diff "$file1" "$file2"; then
echo -e "${name}${blanks:${#name}} OK"
else
echo -e "${name}${blanks:${#name}} FAILED"
Modified: pkg/ProbABEL/checks/test_bt.sh
===================================================================
--- pkg/ProbABEL/checks/test_bt.sh 2014-04-07 15:44:06 UTC (rev 1681)
+++ pkg/ProbABEL/checks/test_bt.sh 2014-04-07 16:34:40 UTC (rev 1682)
@@ -78,10 +78,8 @@
run_diff logist_prob_add.out.txt \
logist_add.out.txt \
- "BT check: prob vs. dose" \
- -I beta_SNP
+ "BT check: prob vs. dose"
run_diff logist_prob_fv_add.out.txt \
logist_fv_add.out.txt \
- "BT check: prob_fv vs. dose_fv" \
- -I beta_SNP
+ "BT check: prob_fv vs. dose_fv"
Modified: pkg/ProbABEL/checks/test_cox.sh
===================================================================
--- pkg/ProbABEL/checks/test_cox.sh 2014-04-07 15:44:06 UTC (rev 1681)
+++ pkg/ProbABEL/checks/test_cox.sh 2014-04-07 16:34:40 UTC (rev 1682)
@@ -59,7 +59,7 @@
>& 3
run_diff coxph_dose_add.out.txt coxph_prob_add.out.txt \
- "pacoxph check: dose vs. prob" -I SNP
+ "pacoxph check: dose vs. prob"
$pacoxph \
Modified: pkg/ProbABEL/checks/test_mms.sh
===================================================================
--- pkg/ProbABEL/checks/test_mms.sh 2014-04-07 15:44:06 UTC (rev 1681)
+++ pkg/ProbABEL/checks/test_mms.sh 2014-04-07 16:34:40 UTC (rev 1682)
@@ -1,7 +1,6 @@
#!/bin/bash
# This script runs checks on ProbABEL's palinear module for
# quantitative traits combined with the mmscore option.
-
echo "Analysis using MMScore..."
scriptdir=$(dirname $0)
@@ -77,10 +76,8 @@
run_diff mmscore_prob_add.out.txt \
mmscore_dose_add.out.txt \
- "mmscore check: prob vs. dose" \
- -I SNP
+ "mmscore check: prob vs. dose"
run_diff mmscore_prob_fv_add.out.txt \
mmscore_dose_fv_add.out.txt \
- "mmscore check: prob_fv vs. dose_fv" \
- -I SNP
+ "mmscore check: prob_fv vs. dose_fv"
Modified: pkg/ProbABEL/checks/test_qt.sh
===================================================================
--- pkg/ProbABEL/checks/test_qt.sh 2014-04-07 15:44:06 UTC (rev 1681)
+++ pkg/ProbABEL/checks/test_qt.sh 2014-04-07 16:34:40 UTC (rev 1682)
@@ -164,7 +164,7 @@
# Remove header from the outputs, because they differ
run_diff linear_base_add.out.txt \
linear_ngp2_add.out.txt \
- "QT check: dose vs. prob (additive model)" -I SNP
+ "QT check: dose vs. prob (additive model)"
for model in add domin over_domin recess 2df; do
run_diff linear_ngp2_${model}.out.txt \
Modified: pkg/ProbABEL/checks/verified_results/height_add.out.txt
===================================================================
--- pkg/ProbABEL/checks/verified_results/height_add.out.txt 2014-04-07 15:44:06 UTC (rev 1681)
+++ pkg/ProbABEL/checks/verified_results/height_add.out.txt 2014-04-07 16:34:40 UTC (rev 1682)
@@ -1,4 +1,4 @@
-name A1 A2 Freq1 MAF Quality Rsq n Mean_predictor_allele chrom position beta_SNP_add sebeta_SNP_add chi2_SNP
+name A1 A2 Freq1 MAF Quality Rsq n Mean_predictor_allele chrom position beta_SNP_addA1 sebeta_SNP_addA1 chi2_SNP_add
rs7247199 GAC A 0.5847 0.415 0.9299 0.8666 182 0.333668 1 204938 20.3847 76.6302 0.0723394
rs8102643 C TGGT 0.5847 0.415 0.9308 0.8685 181 0.808177 1 207859 -25.7632 86.4877 0.0907169
rs8102615 T A 0.5006 0.4702 0.9375 0.8932 182 0.8665 2 211970 -40.7722 102.559 0.161525
Modified: pkg/ProbABEL/checks/verified_results/height_base_add.out.txt
===================================================================
--- pkg/ProbABEL/checks/verified_results/height_base_add.out.txt 2014-04-07 15:44:06 UTC (rev 1681)
+++ pkg/ProbABEL/checks/verified_results/height_base_add.out.txt 2014-04-07 16:34:40 UTC (rev 1682)
@@ -1,4 +1,4 @@
-name A1 A2 Freq1 MAF Quality Rsq n Mean_predictor_allele chrom position beta_SNP_add sebeta_SNP_add chi2_SNP
+name A1 A2 Freq1 MAF Quality Rsq n Mean_predictor_allele chrom position beta_SNP_addA1 sebeta_SNP_addA1 chi2_SNP_add
rs7247199 GAC A 0.5847 0.415 0.9299 0.8666 182 0.333668 19 204938 20.3847 76.6302 0.0723394
rs8102643 C TGGT 0.5847 0.415 0.9308 0.8685 181 0.808177 19 207859 -25.7632 86.4877 0.0907169
rs8102615 T A 0.5006 0.4702 0.9375 0.8932 182 0.8665 19 211970 -40.7722 102.559 0.161525
Modified: pkg/ProbABEL/checks/verified_results/height_ngp2_add.out.txt
===================================================================
--- pkg/ProbABEL/checks/verified_results/height_ngp2_add.out.txt 2014-04-07 15:44:06 UTC (rev 1681)
+++ pkg/ProbABEL/checks/verified_results/height_ngp2_add.out.txt 2014-04-07 16:34:40 UTC (rev 1682)
@@ -1,4 +1,4 @@
-name A1 A2 Freq1 MAF Quality Rsq n Mean_predictor_allele chrom position beta_SNP_addA1 sebeta_SNP_addA1 chi2_SNP_A1
+name A1 A2 Freq1 MAF Quality Rsq n Mean_predictor_allele chrom position beta_SNP_addA1 sebeta_SNP_addA1 chi2_SNP_add
rs7247199 GAC A 0.5847 0.415 0.9299 0.8666 182 0.333668 1 204938 20.3847 76.6302 0.0723394
rs8102643 C TGGT 0.5847 0.415 0.9308 0.8685 181 0.808177 1 207859 -25.7632 86.4877 0.0907169
rs8102615 T A 0.5006 0.4702 0.9375 0.8932 182 0.8665 2 211970 -40.7722 102.559 0.161525
Modified: pkg/ProbABEL/checks/verified_results/height_ngp2_domin.out.txt
===================================================================
--- pkg/ProbABEL/checks/verified_results/height_ngp2_domin.out.txt 2014-04-07 15:44:06 UTC (rev 1681)
+++ pkg/ProbABEL/checks/verified_results/height_ngp2_domin.out.txt 2014-04-07 16:34:40 UTC (rev 1682)
@@ -1,4 +1,4 @@
-name A1 A2 Freq1 MAF Quality Rsq n Mean_predictor_allele chrom position beta_SNP_domA1 sebeta_SNP_domA1 chi2_SNP_domA1
+name A1 A2 Freq1 MAF Quality Rsq n Mean_predictor_allele chrom position beta_SNP_domA1 sebeta_SNP_domA1 chi2_SNP_dom
rs7247199 GAC A 0.5847 0.415 0.9299 0.8666 182 0.333668 1 204938 33.0592 115.333 0.0839902
rs8102643 C TGGT 0.5847 0.415 0.9308 0.8685 181 0.808177 1 207859 -178.139 417.232 0.186313
rs8102615 T A 0.5006 0.4702 0.9375 0.8932 182 0.8665 2 211970 -272.156 640.298 0.18463
Modified: pkg/ProbABEL/checks/verified_results/height_ngp2_over_domin.out.txt
===================================================================
--- pkg/ProbABEL/checks/verified_results/height_ngp2_over_domin.out.txt 2014-04-07 15:44:06 UTC (rev 1681)
+++ pkg/ProbABEL/checks/verified_results/height_ngp2_over_domin.out.txt 2014-04-07 16:34:40 UTC (rev 1682)
@@ -1,4 +1,4 @@
-name A1 A2 Freq1 MAF Quality Rsq n Mean_predictor_allele chrom position beta_SNP_odomA1 sebeta_SNP_odomA1 chi2_SNP_odomA1
+name A1 A2 Freq1 MAF Quality Rsq n Mean_predictor_allele chrom position beta_SNP_odomA1 sebeta_SNP_odomA1 chi2_SNP_odom
rs7247199 GAC A 0.5847 0.415 0.9299 0.8666 182 0.333668 1 204938 80.0802 231.105 0.122726
rs8102643 C TGGT 0.5847 0.415 0.9308 0.8685 181 0.808177 1 207859 28.4969 142.757 0.0407436
rs8102615 T A 0.5006 0.4702 0.9375 0.8932 182 0.8665 2 211970 51.4168 142.042 0.133927
Modified: pkg/ProbABEL/checks/verified_results/height_ngp2_recess.out.txt
===================================================================
--- pkg/ProbABEL/checks/verified_results/height_ngp2_recess.out.txt 2014-04-07 15:44:06 UTC (rev 1681)
+++ pkg/ProbABEL/checks/verified_results/height_ngp2_recess.out.txt 2014-04-07 16:34:40 UTC (rev 1682)
@@ -1,4 +1,4 @@
-name A1 A2 Freq1 MAF Quality Rsq n Mean_predictor_allele chrom position beta_SNP_recA1 sebeta_SNP_recA1 chi2_SNP_recA1
+name A1 A2 Freq1 MAF Quality Rsq n Mean_predictor_allele chrom position beta_SNP_recA1 sebeta_SNP_recA1 chi2_SNP_rec
rs7247199 GAC A 0.5847 0.415 0.9299 0.8666 182 0.333668 1 204938 51.0131 227.445 0.051428
rs8102643 C TGGT 0.5847 0.415 0.9308 0.8685 181 0.808177 1 207859 -28.2737 108.063 0.0699893
rs8102615 T A 0.5006 0.4702 0.9375 0.8932 182 0.8665 2 211970 -45.9705 119.627 0.15093
Modified: pkg/ProbABEL/configure.ac
===================================================================
--- pkg/ProbABEL/configure.ac 2014-04-07 15:44:06 UTC (rev 1681)
+++ pkg/ProbABEL/configure.ac 2014-04-07 16:34:40 UTC (rev 1682)
@@ -21,7 +21,7 @@
AC_PROG_LN_S
if test -z "$CXXFLAGS"; then
# User did not set CXXFLAGS, so we can put in our own defaults
- CXXFLAGS="-g -O2"
+ CXXFLAGS="-g -O2 -DNDEBUG -Wextra"
fi
if test -z "$CPPFLAGS"; then
# User did not set CPPFLAGS, so we can put in our own defaults
Modified: pkg/ProbABEL/examples/mmscore.R
===================================================================
--- pkg/ProbABEL/examples/mmscore.R 2014-04-07 15:44:06 UTC (rev 1681)
+++ pkg/ProbABEL/examples/mmscore.R 2014-04-07 16:34:40 UTC (rev 1682)
@@ -88,85 +88,5 @@
## 2) residuals of the phenotype, which will be the new phenotype that
## ProbABEL will analyse.
-## Mow, go to ProbABEL and start analysis
+## Now, go to ProbABEL and start analysis
-
-
-
-
-##____________________________________________________
-## The following part is for historic purposes only. It is not
-## necessary for using the --mmscore option of ProbABEL's palinear.
-
-## Create test file with genotypes
-
-data_cut <- data[, snps]
-gen_table <- as.numeric(data_cut)
-prob_table <- matrix()
-
-#Replace NA by mean for each snp. NA is forbiden in genotypes in ProbABEL input (!).
-#for(snpnum in 1:dim(gen_table)[2])
-# {
-# mean <- mean(gen_table[,snpnum], na.rm=T)
-# gen_table[is.na(gen_table[,snpnum]),snpnum] <- mean
-# }
-
-
-gen_table_df <- data.frame(gen_table)
-gen_table_df$MLDOSE <- gen_table_df[, 1]
-gen_table_df[,1] <- "MLDOSE"
-colnam <- colnames(gen_table_df)
-
-#colnam[-c(1:length(colnam)-1)] <- colnam[1]
-#colnam[1] <- "MLDOSE"
-
-colnam[length(colnam)] <- colnam[1]
-colnam[1] <- "MLDOSE"
-colnames(gen_table_df) <- colnam
-
-
-rownames <- rownames(gen_table_df)
-rownames <- paste("1->", rownames, sep="")
-rownames(gen_table_df) <- rownames
-
-write.table(file="mmscore_gen.mldose",
- gen_table_df,
- row.names=TRUE,
- col.names=FALSE,
- quote=FALSE,
- na="NaN")
-
-mlinfo <- data.frame(SNP=colnam[2:length(colnam)])
-mlinfo$Al1 <- "A"
-mlinfo$Al2 <- "B"
-mlinfo$Freq1 <- 0.5847
-mlinfo$MAF <- 0.5847
-mlinfo$Quality <- 0.5847
-mlinfo$Rsq <- 0.5847
-
-write.table(mlinfo, "mmscore_gen.mlinfo", row.names=FALSE, quote=FALSE)
-
-## arrange probability-file
-prob_table <- matrix(NA,
- ncol=(dim(gen_table_df)[2]-1) * 2,
- nrow=dim(gen_table_df)[1])
-j <- 1
-for (i in (1:(dim(gen_table_df)[2]-1))) {
- prob_table[,j] <- rep(0, dim(prob_table)[1])
- prob_table[gen_table_df[,i+1]==2, j] <- 1
- prob_table[is.na(gen_table_df[, i+1]), j] <- NA
- j <- j + 1
- prob_table[, j] <- rep(0,dim(prob_table)[1])
- prob_table[gen_table_df[, i+1]==1, j] <- 1
- prob_table[is.na(gen_table_df[, i+1]),j] <- NA
- j <- j + 1
-}
-prob_table_df <- data.frame(MLPROB="MLPROB", prob_table)
-rownames(prob_table_df) <- rownames(gen_table_df)
-
-write.table(file="mmscore_gen.mlprob",
- prob_table_df,
- row.names=TRUE,
- col.names=FILE,
- quote=FILE,
- na="NaN")
Modified: pkg/ProbABEL/src/cholesky.cpp
===================================================================
--- pkg/ProbABEL/src/cholesky.cpp 2014-04-07 15:44:06 UTC (rev 1681)
+++ pkg/ProbABEL/src/cholesky.cpp 2014-04-07 16:34:40 UTC (rev 1682)
@@ -100,8 +100,8 @@
}
int n = matrix.ncol;
- register double temp;
- register int i, j, k;
+ double temp;
+ int i, j, k;
/*
** invert the cholesky in the lower triangle
Modified: pkg/ProbABEL/src/command_line_settings.cpp
===================================================================
--- pkg/ProbABEL/src/command_line_settings.cpp 2014-04-07 15:44:06 UTC (rev 1681)
+++ pkg/ProbABEL/src/command_line_settings.cpp 2014-04-07 16:34:40 UTC (rev 1682)
@@ -120,12 +120,12 @@
{
return outfilename;
}
+//TODO(unknown) This function is not used. Remove in near future
+//char* cmdvars::getProgramName() const
+//{
+// return program_name;
+//}
-char* cmdvars::getProgramName() const
-{
- return program_name;
-}
-
int cmdvars::getRobust() const
{
return robust;
Modified: pkg/ProbABEL/src/command_line_settings.h
===================================================================
--- pkg/ProbABEL/src/command_line_settings.h 2014-04-07 15:44:06 UTC (rev 1681)
+++ pkg/ProbABEL/src/command_line_settings.h 2014-04-07 16:34:40 UTC (rev 1682)
@@ -116,7 +116,8 @@
int getNoutcomes() const;
int getNpeople() const;
string getOutfilename() const;
- char* getProgramName() const;
+//TODO(unknown) This function is not used. Remove in near future
+// char* getProgramName() const;
int getRobust() const;
int getScore() const;
string getSep() const;
Modified: pkg/ProbABEL/src/coxph_data.cpp
===================================================================
--- pkg/ProbABEL/src/coxph_data.cpp 2014-04-07 15:44:06 UTC (rev 1681)
+++ pkg/ProbABEL/src/coxph_data.cpp 2014-04-07 16:34:40 UTC (rev 1682)
@@ -81,10 +81,7 @@
freq = 0;
masked_data = new unsigned short int[nids];
- for (int i = 0; i < nids; i++)
- {
- masked_data[i] = obj.masked_data[i];
- }
+ std::copy(obj.masked_data, obj.masked_data+nids,masked_data);
}
coxph_data::coxph_data(phedata &phed, gendata &gend, const int snpnum)
@@ -94,11 +91,9 @@
nids = gend.nids;
masked_data = new unsigned short int[nids];
- for (int i = 0; i < nids; i++)
- {
- masked_data[i] = 0;
- }
+ std::fill(masked_data,masked_data+nids,0);
+
ngpreds = gend.ngpreds;
if (snpnum >= 0)
{
@@ -167,11 +162,12 @@
// sort by time
double *tmptime = new double[nids];
int *passed_sorted = new int[nids];
+ std::fill (passed_sorted,passed_sorted+nids,0);
+
for (int i = 0; i < nids; i++)
{
tmptime[i] = stime[i];
- passed_sorted[i] = 0;
}
qsort(tmptime, nids, sizeof(double), cmpfun);
@@ -239,9 +235,7 @@
for (int j = 0; j < ngpreds; j++) {
double *snpdata = new double[nids];
- for (int i = 0; i < nids; i++) {
- masked_data[i] = 0;
- }
+ std::fill (masked_data,masked_data+nids,0);
gend->get_var(snpnum * ngpreds + j, snpdata);
@@ -316,15 +310,9 @@
coxph_data to; // = coxph_data(*this);
// filter missing data
- int nmeasured = 0;
- for (int i = 0; i < nids; i++)
- {
- if (masked_data[i] == 0)
- {
- nmeasured++;
- }
- }
+ int nmeasured=std::count (masked_data, masked_data+nids, 0);
+
to.nids = nmeasured;
to.ncov = ncov;
to.ngpreds = ngpreds;
@@ -358,11 +346,7 @@
//delete [] to.masked_data;
to.masked_data = new unsigned short int[to.nids];
- for (int i = 0; i < to.nids; i++)
- {
- to.masked_data[i] = 0;
- }
-
+ std::fill(to.masked_data,to.masked_data+to.nids,0);
return (to);
}
@@ -379,7 +363,7 @@
}
-void coxph_reg::estimate(coxph_data &cdatain, const int verbose,
+void coxph_reg::estimate(coxph_data &cdatain,
int maxiter, double eps,
double tol_chol, const int model,
const int interaction, const int ngpreds,
Modified: pkg/ProbABEL/src/coxph_data.h
===================================================================
--- pkg/ProbABEL/src/coxph_data.h 2014-04-07 15:44:06 UTC (rev 1681)
+++ pkg/ProbABEL/src/coxph_data.h 2014-04-07 16:34:40 UTC (rev 1682)
@@ -89,7 +89,7 @@
int niter;
coxph_reg(coxph_data &cdatain);
- void estimate(coxph_data &cdatain, const int verbose, int maxiter,
+ void estimate(coxph_data &cdatain, int maxiter,
double eps, double tol_chol, const int model,
const int interaction, const int ngpreds, const bool iscox,
const int nullmodel, const mlinfo &snpinfo, const int cursnp);
Modified: pkg/ProbABEL/src/data.cpp
===================================================================
--- pkg/ProbABEL/src/data.cpp 2014-04-07 15:44:06 UTC (rev 1681)
+++ pkg/ProbABEL/src/data.cpp 2014-04-07 16:34:40 UTC (rev 1682)
@@ -47,49 +47,45 @@
#endif
#include "utilities.h"
+//TODO(unknown) This function is not used. Remove in near future
+//unsigned int Nmeasured(char * fname, int nphenocols, int npeople)
+//{
+//// first pass -- find unmeasured people
+// std::ifstream infile(fname);
+// if (!infile)
+// {
+// std::cerr << "Nmeasured: cannot open file " << fname << endl;
+// }
+// char tmp[100];
+//
+// for (int i = 0; i < nphenocols; i++)
+// {
+// infile >> tmp;
+// }
+//
+// unsigned short int * allmeasured = new unsigned short int[npeople];
+// int nids = 0;
+// for (int i = 0; i < npeople; i++)
+// {
+// allmeasured[i] = 1;
+// infile >> tmp;
+// for (int j = 1; j < nphenocols; j++)
+// {
+// infile >> tmp;
+// if (tmp[0] == 'N' || tmp[0] == 'n')
+// allmeasured[i] = 0;
+// }
+// if (allmeasured[i] == 1)
+// nids++;
+// }
+// infile.close();
+//
+// delete[] allmeasured;
+//
+// return (nids);
+//}
-unsigned int Nmeasured(char * fname, int nphenocols, int npeople)
-{
-//TODO: unused variables remove them for good if there is no reason to keep them
-//int ncov = nphenocols - 2;
-//int nids_all = npeople;
-// first pass -- find unmeasured people
- std::ifstream infile(fname);
- if (!infile)
- {
- std::cerr << "Nmeasured: cannot open file " << fname << endl;
- }
- char tmp[100];
-
- for (int i = 0; i < nphenocols; i++)
- {
- infile >> tmp;
- }
-
- unsigned short int * allmeasured = new unsigned short int[npeople];
- int nids = 0;
- for (int i = 0; i < npeople; i++)
- {
- allmeasured[i] = 1;
- infile >> tmp;
- for (int j = 1; j < nphenocols; j++)
- {
- infile >> tmp;
- if (tmp[0] == 'N' || tmp[0] == 'n')
- allmeasured[i] = 0;
- }
- if (allmeasured[i] == 1)
- nids++;
- }
- infile.close();
-
- delete[] allmeasured;
-
- return (nids);
-}
-
-
/**
* Read SNP information from an mlinfo file generated by the
* imputation software.
@@ -269,5 +265,3 @@
{
return matrix;
}
-
-//________________________________________Maksim_end
Modified: pkg/ProbABEL/src/data.h
===================================================================
--- pkg/ProbABEL/src/data.h 2014-04-07 15:44:06 UTC (rev 1681)
+++ pkg/ProbABEL/src/data.h 2014-04-07 16:34:40 UTC (rev 1682)
@@ -31,8 +31,8 @@
#include <string>
extern bool is_interaction_excluded;
-
-unsigned int Nmeasured(char * fname, int nphenocols, int npeople);
+//TODO(unknown) This function is not used. Remove in near future
+//unsigned int Nmeasured(char * fname, int nphenocols, int npeople);
#include "phedata.h"
#include "gendata.h"
@@ -80,4 +80,4 @@
~InvSigma();
};
-#endif /* DATA_H_ */
+#endif//DATA_H_
Modified: pkg/ProbABEL/src/fvlib
===================================================================
--- pkg/ProbABEL/src/fvlib 2014-04-07 15:44:06 UTC (rev 1681)
+++ pkg/ProbABEL/src/fvlib 2014-04-07 16:34:40 UTC (rev 1682)
@@ -1 +1 @@
-link ../../filevector/fvlib
\ No newline at end of file
+link include/filevector/fvlib
\ No newline at end of file
Modified: pkg/ProbABEL/src/gendata.cpp
===================================================================
--- pkg/ProbABEL/src/gendata.cpp 2014-04-07 15:44:06 UTC (rev 1681)
+++ pkg/ProbABEL/src/gendata.cpp 2014-04-07 16:34:40 UTC (rev 1682)
@@ -28,6 +28,7 @@
#include <string>
#include <errno.h>
+#include <limits>
#include "gendata.h"
#include "fvlib/FileVector.h"
#if EIGEN
@@ -39,6 +40,69 @@
#endif
#include "utilities.h"
+
+void gendata::mldose_line_to_matrix(int k,const char *all_numbers,int amount_of_numbers){
+ int j = 0;
+ //check if not a null pointer
+ if (!*all_numbers){
+ perror("Error while reading genetic data (expected pointer to char but found a null pointer)");
+ exit(EXIT_FAILURE);
+ }
+ while (j<amount_of_numbers)
+ {
+ double result = 0;
+ //skip whitespace
+ while (*all_numbers == ' ')
+ {
+ all_numbers++;
+ }
+ //check NaN (right now checks only first character)
+ //TODO: make catching of NaN more rigid
+ if (*all_numbers == 'N')
+ {
+ result = std::numeric_limits<double>::quiet_NaN();
+ //skip other characters of NaN
+ while ((*all_numbers == 'a') | (*all_numbers == 'N'))
+ {
+ all_numbers++;
+ }
+ }
+ else
+ {
+ int sign = 0;
+ //set sign to -1 if negative: multiply by sign just before return
+ if (*all_numbers == '-')
+ {
+ all_numbers++;
+ sign = -1;
+ }
+ //read digits before dot
+ while (*all_numbers <= '9' && *all_numbers >= '0')
+ {
+ result = result * 10 + (*all_numbers++ - '0');
+ }
+ //read digit after dot
+ if (*all_numbers == '.')
+ {
+ double decimal_counter = 1.0;
+ all_numbers++;
+ while (*all_numbers <= '9' && *all_numbers >= '0')
+ {
+ decimal_counter *= 0.1;
+ result += (*all_numbers++ - '0') * decimal_counter;
+ }
+ }
+ //correct for negative number
+ if (sign == -1)
+ {
+ result = sign * result;
+ }
+ }
+ G.put(result, k, j);
+ j++;
+ }
+}
+
void gendata::get_var(int var, double * data)
{
// Read the genetic data for SNP 'var' and store in the array 'data'
@@ -165,17 +229,19 @@
DAG = NULL;
// int nids_all = npeople;
+
G.reinit(nids, (nsnps * ngpreds));
std::ifstream infile;
+ infile.open(fname);
- infile.open(fname);
if (!infile)
{
std::cerr << "gendata: cannot open file " << fname << endl;
}
std::string tmpid, tmpstr;
+ char inStr[8];
int k = 0;
for (unsigned int i = 0; i < npeople; i++)
@@ -191,7 +257,7 @@
size_t strpos = tmpstr.find("->");
if (strpos != string::npos)
{
- tmpid = tmpstr.substr(strpos+2, string::npos);
+ tmpid = tmpstr.substr(strpos + 2, string::npos);
}
else
{
@@ -200,8 +266,8 @@
if (tmpid != idnames[k])
{
cerr << "phenotype file and dose or probability file "
- << "did not match at line " << i + 2 << " (" << tmpid
- << " != " << idnames[k] << ")" << endl;
+ << "did not match at line " << i + 2 << " ("
+ << tmpid << " != " << idnames[k] << ")" << endl;
infile.close();
exit(1);
}
@@ -212,46 +278,58 @@
infile >> tmpstr;
}
- for (unsigned int j = 0; j < (nsnps * ngpreds); j++)
+ int oldstyle = 0;
+ if (oldstyle == 1)
{
- if (infile.good())
+ for (unsigned int j = 0; j < (nsnps * ngpreds); j++)
{
- infile >> tmpstr;
- // tmpstr contains the dosage/probability in
- // string form. Convert it to double (if tmpstr is
- // NaN it will be set to nan).
- double dosage;
- char *endptr;
- errno = 0; // To distinguish success/failure
- // after strtod()
- dosage = strtod(tmpstr.c_str(), &endptr);
+ if (infile.good())
+ {
+ infile >> inStr;
+ // tmpstr contains the dosage/probability in
+ // string form. Convert it to double (if tmpstr is
+ // NaN it will be set to nan).
+ double dosage;
+ char *endptr;
+ errno = 0; // To distinguish success/failure
+ // after strtod()
- if ((errno == ERANGE &&
- (dosage == HUGE_VALF || dosage == HUGE_VALL))
- || (errno != 0 && dosage == 0)) {
- perror("Error while reading genetic data (strtod)");
- exit(EXIT_FAILURE);
- }
+ dosage = strtod(inStr, &endptr);
+ if ((errno == ERANGE
+ && (dosage == HUGE_VALF || dosage == HUGE_VALL))
+ || (errno != 0 && dosage == 0))
+ {
+ perror("Error while reading genetic data (strtod)");
+ exit(EXIT_FAILURE);
+ }
- if (endptr == tmpstr.c_str()) {
- cerr << "No digits were found while reading genetic data"
- << " (individual " << i + 1
- << ", position " << j + 1 << ")"
- << endl;
- exit(EXIT_FAILURE);
+ if (endptr == tmpstr.c_str())
+ {
+ cerr
+ << "No digits were found while reading genetic data"
+ << " (individual " << i + 1 << ", position "
+ << j + 1 << ")" << endl;
+ exit(EXIT_FAILURE);
+ }
+ /* If we got here, strtod() successfully parsed a number */
+ G.put(dosage, k, j);
}
-
- /* If we got here, strtod() successfully parsed a number */
- G.put(dosage, k, j);
+ else
+ {
+ std::cerr << "cannot read dose-file: " << fname
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/genabel -r 1682
More information about the Genabel-commits
mailing list