[Genabel-commits] r1231 - in pkg/ProbABEL: . doc examples src tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon May 20 01:05:40 CEST 2013
Author: lckarssen
Date: 2013-05-20 01:05:39 +0200 (Mon, 20 May 2013)
New Revision: 1231
Added:
pkg/ProbABEL/tests/test_cox_ngpreds2/
Modified:
pkg/ProbABEL/configure.ac
pkg/ProbABEL/doc/pacoxph.1
pkg/ProbABEL/examples/Makefile.am
pkg/ProbABEL/examples/example_cox.sh
pkg/ProbABEL/examples/example_mms.sh
pkg/ProbABEL/src/coxfit2.c
pkg/ProbABEL/src/coxph_data.cpp
pkg/ProbABEL/src/coxph_data.h
pkg/ProbABEL/src/eigen_mematrix.cpp
pkg/ProbABEL/src/main.cpp
pkg/ProbABEL/src/reg1.cpp
Log:
ProbABEL trunk: Imported changes from the PAcoxfix branch. This change finally makes the Cox PH model work again (it was not working (correctly)) since the implementation of the filevector lib.
This should fix longstanding bug #1130.
These changes were partially sponsored by the Erasmus Medical Center, Rotterdam (the group of prof. C.M. van Duijn) and the Radboud Medical Center, Nijmegen (the group of prof. Kiemeney); as well as a considerable amount of my own time.
Here is a break-down of the changes per file.
examples/example_cox.sh: Bring the example script for Cox up to date with the others.
examples/Makefile.am: Enable the Cox tests and list the files that are created by the example script.
examples/example_mms.sh: Unrelated change, add "dose" to mmscore output files made with dose input.
tests/test_cox_ngpreds2: Add R test scripts to validate the output of the linear and Cox modules. These tests are not run automatically by "make check".
configure.ac: Now building pacoxph is enabled by default.
doc/pacoxph.1: Remove warning message about pacoxph being buggy from the man page.
src/reg1.cpp: An important part of this fix is in lines 209-210, where an argument was missing. As a result the 'iscox' variable wasn't passed in the right slot. This went unnoticed because the last variable ("nullmodel") had a default value.
src/coxph_data.h: only cosmetic changes to the code layout.
src/coxph_data.cpp: many code layout improvements. Important part of the fix in lines 257-262.
src/eigen_mematrix.cpp: Important part of the fix in line 316, reordering the data went wrong here.
src/coxfit2.c: Fixing a few (unrelated to the real Cox problems) memory leaks in lines 302-305, 420-423 and 461-464.
src/main.cpp: A bunch of code layout improvements. The important part of the fix is in lines 613-616 where the "iscox" argument (set to true) wasn't passed on correctly, just like in reg1.cpp.
Modified: pkg/ProbABEL/configure.ac
===================================================================
--- pkg/ProbABEL/configure.ac 2013-05-19 21:27:46 UTC (rev 1230)
+++ pkg/ProbABEL/configure.ac 2013-05-19 23:05:39 UTC (rev 1231)
@@ -111,10 +111,9 @@
# Since pacoxph is buggy it needs to be enabled explixitly
AC_ARG_ENABLE([pacoxph],
- [AS_HELP_STRING([--enable-pacoxph], [build the pacoxph binary (still \
-contains lots of bugs)])],
- [pacoxph=${enableval}],
- [pacoxph=no])
+ [AS_HELP_STRING([--disable-pacoxph], [disable building the pacoxph program])],
+ [pacoxph=no],
+ [pacoxph=yes])
if test "x$pacoxph" = "xyes"; then
AC_MSG_NOTICE([building of pacoxph is enabled])
Modified: pkg/ProbABEL/doc/pacoxph.1
===================================================================
--- pkg/ProbABEL/doc/pacoxph.1 2013-05-19 21:27:46 UTC (rev 1230)
+++ pkg/ProbABEL/doc/pacoxph.1 2013-05-19 23:05:39 UTC (rev 1231)
@@ -1,4 +1,4 @@
-.TH pacoxph 1 "23 February 2012"
+.TH pacoxph 1 "08 February 2013"
.SH NAME
pacoxph \- Perform Genome-Wide Association Analysis using a linear model
.SH SYNOPSIS
@@ -78,8 +78,6 @@
.SH "SEE ALSO"
palinear(1), palogist(1)
.SH BUGS
-Unfortunately
-.B pacoxph
-is in a buggy state at the moment. It cannot use files in DatABEL format.
+None known at the moment. The bugtracker is located at https://r-forge.r-project.org/tracker/?group_id=505
.SH AUTHORS
Lennart C. Karssen
Modified: pkg/ProbABEL/examples/Makefile.am
===================================================================
--- pkg/ProbABEL/examples/Makefile.am 2013-05-19 21:27:46 UTC (rev 1230)
+++ pkg/ProbABEL/examples/Makefile.am 2013-05-19 23:05:39 UTC (rev 1231)
@@ -46,7 +46,9 @@
TESTS_ENVIRONMENT = bash
check_SCRIPTS = example_bt.sh example_qt.sh example_mms.sh
-# example_cox.sh: don't run this test. It takes several hours...
+if BUILD_pacoxph
+check_SCRIPTS += example_cox.sh
+endif
TESTS = $(check_SCRIPTS)
@@ -98,13 +100,18 @@
height_ngp2_allcov_fv_2df.out.txt height_ngp2_robust_fv_add.out.txt \
height_ngp2_allcov_2df.out.txt
-cleanfiles_mms = mmscore_add.out.txt mmscore_fv_add.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 \
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_add.out.txt coxph_fv_add.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 = $(cleanfiles_bt) $(cleanfiles_qt) $(cleanfiles_mms) $(cleanfiles_cox)
Modified: pkg/ProbABEL/examples/example_cox.sh
===================================================================
--- pkg/ProbABEL/examples/example_cox.sh 2013-05-19 21:27:46 UTC (rev 1230)
+++ pkg/ProbABEL/examples/example_cox.sh 2013-05-19 23:05:39 UTC (rev 1231)
@@ -1,36 +1,69 @@
-echo "analysing Cox model"
+#!/bin/sh
+# This script runs checks on ProbABEL's pacoxph module
+
+echo "Analysing Cox model..."
if [ -z ${srcdir} ]; then
srcdir="."
fi
+. ${srcdir}/run_diff.sh
+
+# 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
+
../src/pacoxph \
-p ${srcdir}/coxph_data.txt \
-d ${srcdir}/test.mldose \
-i ${srcdir}/test.mlinfo \
-m ${srcdir}/test.map \
-c 19 \
- -o coxph
+ -o coxph_dose \
+ >& 3
-# ../src/pacoxph \
-# -p ${srcdir}/coxph_data.txt \
-# -d ${srcdir}/test.dose.fvi \
-# -i ${srcdir}/test.mlinfo \
-# -m ${srcdir}/test.map \
-# -c 19 \
-# -o coxph_fv
+../src/pacoxph \
+ -p ${srcdir}/coxph_data.txt \
+ -d ${srcdir}/test.dose.fvi \
+ -i ${srcdir}/test.mlinfo \
+ -m ${srcdir}/test.map \
+ -c 19 \
+ -o coxph_dose_fv \
+ >& 3
-# echo "pacoxph check: dose vs. dose_fv"
-# diff coxph_add.out.txt coxph_fv_add.out.txt
+run_diff coxph_dose_add.out.txt coxph_dose_fv_add.out.txt \
+ "pacoxph check: dose vs. dose_fv"
-# ../src/pacoxph \
-# -p ${srcdir}/coxph_data.txt \
-# -d ${srcdir}/test.mlprob \
-# -i ${srcdir}/test.mlinfo \
-# -m ${srcdir}/test.map \
-# --ngpreds=2 \
-# -c 19 \
-# -o coxph_prob
+../src/pacoxph \
+ -p ${srcdir}/coxph_data.txt \
+ -d ${srcdir}/test.mlprob \
+ -i ${srcdir}/test.mlinfo \
+ -m ${srcdir}/test.map \
+ --ngpreds=2 \
+ -c 19 \
+ -o coxph_prob \
+ >& 3
-# echo "pacoxph check: dose vs. prob"
-# diff coxph_add.out.txt coxph_prob_add.out.txt
+run_diff coxph_dose_add.out.txt coxph_prob_add.out.txt \
+ "pacoxph check: dose vs. prob" -I SNP
+
+
+../src/pacoxph \
+ -p ${srcdir}/coxph_data.txt \
+ -d ${srcdir}/test.prob.fvi \
+ -i ${srcdir}/test.mlinfo \
+ -m ${srcdir}/test.map \
+ --ngpreds=2 \
+ -c 19 \
+ -o coxph_prob_fv \
+ >& 3
+
+for model in add domin recess over_domin 2df; do
+ run_diff coxph_prob_${model}.out.txt \
+ coxph_prob_fv_${model}.out.txt \
+ "pacoxph check ($model model): prob vs. prob_fv"
+done
Modified: pkg/ProbABEL/examples/example_mms.sh
===================================================================
--- pkg/ProbABEL/examples/example_mms.sh 2013-05-19 21:27:46 UTC (rev 1230)
+++ pkg/ProbABEL/examples/example_mms.sh 2013-05-19 23:05:39 UTC (rev 1231)
@@ -21,7 +21,7 @@
-i ${srcdir}/mmscore_gen.mlinfo \
-d ${srcdir}/mmscore_gen.mldose \
--sep="," \
- -o mmscore \
+ -o mmscore_dose \
--mmscore ${srcdir}/mmscore_InvSigma_aj.sex.age.dat \
>& 3
@@ -30,13 +30,13 @@
-i ${srcdir}/mmscore_gen.mlinfo \
-d ${srcdir}/mmscore_gen.dose.fvi \
--sep="," \
- -o mmscore_fv \
+ -o mmscore_dose_fv \
--mmscore ${srcdir}/mmscore_InvSigma_aj.sex.age.dat \
>& 3
-run_diff mmscore_add.out.txt \
- mmscore_fv_add.out.txt \
+run_diff mmscore_dose_add.out.txt \
+ mmscore_dose_fv_add.out.txt \
"mmscore check: dose vs. dose_fv"
Modified: pkg/ProbABEL/src/coxfit2.c
===================================================================
--- pkg/ProbABEL/src/coxfit2.c 2013-05-19 21:27:46 UTC (rev 1230)
+++ pkg/ProbABEL/src/coxfit2.c 2013-05-19 23:05:39 UTC (rev 1231)
@@ -181,7 +181,7 @@
a2[i] = 0;
for (j = 0; j < nvar; j++)
{
-// std::cout << i << " " << j << "\n";
+// std::cout << i << " " << j << "\n";
cmat[i][j] = 0;
cmat2[i][j] = 0;
}
@@ -194,7 +194,7 @@
for (i = 0; i < nvar; i++)
{
zbeta += beta[i] * covar[i][person];
-// std::cout << zbeta << " " << beta[i] << " " << covar[i][person] << "\n";
+// std::cout << zbeta << " " << beta[i] << " " << covar[i][person] << "\n";
}
risk = exp(zbeta) * weights[person];
@@ -242,7 +242,7 @@
{
temp2 = (a[i] - temp * a2[i]) / d2;
u[i] -= wtave[person] * temp2;
-// std::cout << i << " " << wtave[person] << " " << temp2 << " " << a[i] << " " << a2[i] << "\n";
+// std::cout << i << " " << wtave[person] << " " << temp2 << " " << a[i] << " " << a2[i] << "\n";
for (j = 0; j <= i; j++)
imat[j][i] += wtave[person]
* ((cmat[i][j] - temp * cmat2[i][j]) / d2
@@ -250,7 +250,7 @@
}
}
// std::cout << person << "\n";
-// print_prematrix(u,1,nvar);
+// print_prematrix(u,1,nvar);
efron_wt = 0;
for (i = 0; i < nvar; i++)
{
@@ -298,6 +298,11 @@
for (i = 1; i < nvar; i++)
for (j = 0; j < i; j++)
imat[i][j] = imat[j][i];
+
+ free(covar);
+ free(imat);
+ free(cmat);
+ free(cmat2);
return; /* and we leave the old beta in peace */
}
@@ -411,6 +416,11 @@
if (halving == 1)
*flag = 1000; /*didn't converge after all */
*maxiter = iter;
+
+ free(covar);
+ free(imat);
+ free(cmat);
+ free(cmat2);
return;
}
@@ -446,5 +456,11 @@
for (i = 0; i < nvar; i++)
beta[i] = newbeta[i];
*flag = 1000;
+
+
+ free(covar);
+ free(imat);
+ free(cmat);
+ free(cmat2);
return;
}
Modified: pkg/ProbABEL/src/coxph_data.cpp
===================================================================
--- pkg/ProbABEL/src/coxph_data.cpp 2013-05-19 21:27:46 UTC (rev 1230)
+++ pkg/ProbABEL/src/coxph_data.cpp 2013-05-19 23:05:39 UTC (rev 1231)
@@ -46,18 +46,21 @@
coxph_data::coxph_data(const coxph_data &obj) : sstat(obj.sstat)
{
- nids = obj.nids;
- ncov = obj.ncov;
- ngpreds = obj.ngpreds;
- weights = obj.weights;
- stime = obj.stime;
- offset = obj.offset;
- strata = obj.strata;
- X = obj.X;
- order = obj.order;
+ nids = obj.nids;
+ ncov = obj.ncov;
+ ngpreds = obj.ngpreds;
+ weights = obj.weights;
+ stime = obj.stime;
+ offset = obj.offset;
+ strata = obj.strata;
+ X = obj.X;
+ order = obj.order;
masked_data = new unsigned short int[nids];
+
for (int i = 0; i < nids; i++)
+ {
masked_data[i] = 0;
+ }
}
coxph_data::coxph_data(phedata &phed, gendata &gend, int snpnum)
@@ -80,7 +83,6 @@
exit(1);
}
- // X.reinit(nids,(ncov+1));
X.reinit(nids, ncov);
stime.reinit(nids, 1);
sstat.reinit(nids, 1);
@@ -103,10 +105,13 @@
}
for (int j = 0; j < phed.ncov; j++)
+ {
for (int i = 0; i < nids; i++)
X.put((phed.X).get(i, j), i, j);
+ }
if (snpnum > 0)
+ {
for (int j = 0; j < ngpreds; j++)
{
double snpdata[nids];
@@ -114,17 +119,15 @@
for (int i = 0; i < nids; i++)
X.put(snpdata[i], i, (ncov - ngpreds + j));
}
+ }
- // for (int i=0;i<nids;i++)
- // for (int j=0;j<ngpreds;j++)
- // X.put((gend.G).get(i,(snpnum*ngpreds+j)),i,(ncov-ngpreds+j));
-
for (int i = 0; i < nids; i++)
{
weights[i] = 1.0;
offset[i] = 0.0;
strata[i] = 0;
}
+
// sort by time
double tmptime[nids];
int passed_sorted[nids];
@@ -134,12 +137,16 @@
tmptime[i] = stime[i];
passed_sorted[i] = 0;
}
+
qsort(tmptime, nids, sizeof(double), cmpfun);
+
for (int i = 0; i < nids; i++)
{
int passed = 0;
for (int j = 0; j < nids; j++)
+ {
if (tmptime[j] == stime[i])
+ {
if (!passed_sorted[j])
{
order[i] = j;
@@ -147,25 +154,29 @@
passed = 1;
break;
}
+ }
+ }
if (passed != 1)
{
std::cerr << "cannot recover element " << i << "\n";
exit(1);
}
}
- stime = reorder(stime, order);
- sstat = reorder(sstat, order);
+ stime = reorder(stime, order);
+ sstat = reorder(sstat, order);
weights = reorder(weights, order);
- strata = reorder(strata, order);
- offset = reorder(offset, order);
- X = reorder(X, order);
- X = transpose(X);
- // X.print();
- // offset.print();
- // weights.print();
- // stime.print();
- // sstat.print();
+ strata = reorder(strata, order);
+ offset = reorder(offset, order);
+ X = reorder(X, order);
+ X = transpose(X);
+
+ // X.print();
+ // offset.print();
+ // weights.print();
+ // stime.print();
+ // sstat.print();
}
+
void coxph_data::update_snp(gendata &gend, int snpnum)
{
/**
@@ -186,7 +197,9 @@
{
double snpdata[nids];
for (int i = 0; i < nids; i++)
+ {
masked_data[i] = 0;
+ }
gend.get_var(snpnum * ngpreds + j, snpdata);
@@ -213,22 +226,21 @@
coxph_data coxph_data::get_unmasked_data()
{
-// std::cout << " !!! in get_unmasked_data !!! ";
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++;
+ }
+ }
to.nids = nmeasured;
-// std::cout << "nmeasured=" << nmeasured << "\n";
to.ncov = ncov;
-// std::cout << "ncov=" << ncov << "\n";
to.ngpreds = ngpreds;
int dim1X = X.nrow;
-// std::cout << "X.ncol=" << X.ncol << "\n";
-// std::cout << "X.nrow=" << X.nrow << "\n";
(to.weights).reinit(to.nids, 1);
(to.stime).reinit(to.nids, 1);
(to.sstat).reinit(to.nids, 1);
@@ -237,34 +249,32 @@
(to.order).reinit(to.nids, 1);
(to.X).reinit(dim1X, to.nids);
-// std::cout << "(to.X).ncol=" << (to.X).ncol << "\n";
-// std::cout << "(to.X).nrow=" << (to.X).nrow << "\n";
-// std::cout << " !!! just before cycle !!! ";
int j = 0;
for (int i = 0; i < nids; i++)
{
-// std::cout << nids << " " << i << " " << masked_data[i] << "\n";
if (masked_data[i] == 0)
{
- (to.weights).put(weights.get(i, 1), j, 1);
- (to.stime).put(stime.get(i, 1), j, 1);
- (to.sstat).put(sstat.get(i, 1), j, 1);
- (to.offset).put(offset.get(i, 1), j, 1);
- (to.strata).put(strata.get(i, 1), j, 1);
- (to.order).put(order.get(i, 1), j, 1);
+ (to.weights).put(weights.get(i, 0), j, 0);
+ (to.stime).put(stime.get(i, 0), j, 0);
+ (to.sstat).put(sstat.get(i, 0), j, 0);
+ (to.offset).put(offset.get(i, 0), j, 0);
+ (to.strata).put(strata.get(i, 0), j, 0);
+ (to.order).put(order.get(i, 0), j, 0);
for (int nc = 0; nc < dim1X; nc++)
+ {
(to.X).put(X.get(nc, i), nc, j);
+ }
j++;
}
}
-// std::cout << " !!! just after cycle !!! ";
//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::cout << "get_unmasked: " << to.nids << " "
- // << dim2X << " " << dim2Y << "\n";
+ }
+
return (to);
}
@@ -287,7 +297,7 @@
coxph_data cdata = cdatain.get_unmasked_data();
mematrix<double> X = t_apply_model(cdata.X, model, interaction, ngpreds,
iscox, nullmodel);
- // X.print();
+
int length_beta = X.nrow;
beta.reinit(length_beta, 1);
sebeta.reinit(length_beta, 1);
@@ -307,19 +317,29 @@
int flag;
double sctest = 1.0;
- //TODO(maarten): this function works only in combination with eigen
- // remove .data() from arguments to make is available under old matrix
+ // When using Eigen coxfit2 needs to be called in a slightly
+ // different way (i.e. the .data()-part needs to be added).
+#if EIGEN
coxfit2(&maxiter, &cdata.nids, &X.nrow, cdata.stime.data.data(),
cdata.sstat.data.data(), X.data.data(), newoffset.data.data(),
cdata.weights.data.data(), cdata.strata.data.data(),
means.data.data(), beta.data.data(), u.data.data(),
imat.data.data(), loglik_int, &flag, work, &eps, &tol_chol,
&sctest);
+#else
+ coxfit2(&maxiter, &cdata.nids, &X.nrow, cdata.stime.data,
+ cdata.sstat.data, X.data, newoffset.data,
+ cdata.weights.data, cdata.strata.data,
+ means.data, beta.data, u.data,
+ imat.data, loglik_int, &flag, work, &eps, &tol_chol,
+ &sctest);
+#endif
for (int i = 0; i < X.nrow; i++)
{
sebeta[i] = sqrt(imat.get(i, i));
}
+
loglik = loglik_int[1];
niter = maxiter;
}
Modified: pkg/ProbABEL/src/coxph_data.h
===================================================================
--- pkg/ProbABEL/src/coxph_data.h 2013-05-19 21:27:46 UTC (rev 1230)
+++ pkg/ProbABEL/src/coxph_data.h 2013-05-19 23:05:39 UTC (rev 1231)
@@ -23,6 +23,7 @@
class coxph_data {
public:
coxph_data get_unmasked_data();
+
coxph_data()
{
nids = 0;
@@ -30,6 +31,7 @@
ngpreds = 0;
masked_data = NULL;
}
+
coxph_data(const coxph_data &obj);
coxph_data(phedata &phed, gendata &gend, int snpnum);
void update_snp(gendata &gend, int snpnum);
@@ -40,11 +42,11 @@
int ngpreds;
mematrix<double> weights;
mematrix<double> stime;
- mematrix<int> sstat;
+ mematrix<int> sstat;
mematrix<double> offset;
- mematrix<int> strata;
+ mematrix<int> strata;
mematrix<double> X;
- mematrix<int> order;
+ mematrix<int> order;
unsigned short int * masked_data;
};
Modified: pkg/ProbABEL/src/eigen_mematrix.cpp
===================================================================
--- pkg/ProbABEL/src/eigen_mematrix.cpp 2013-05-19 21:27:46 UTC (rev 1230)
+++ pkg/ProbABEL/src/eigen_mematrix.cpp 2013-05-19 23:05:39 UTC (rev 1231)
@@ -313,7 +313,7 @@
for (int i = 0; i < temp.nrow; i++)
{
source = order.data(i, 0);
- temp.data.row(i) = M.data.row(source);
+ temp.data.row(source) = M.data.row(i);
}
return temp;
}
Modified: pkg/ProbABEL/src/main.cpp
===================================================================
--- pkg/ProbABEL/src/main.cpp 2013-05-19 21:27:46 UTC (rev 1230)
+++ pkg/ProbABEL/src/main.cpp 2013-05-19 23:05:39 UTC (rev 1231)
@@ -59,18 +59,18 @@
if (csnp == 0)
{
std::cout << "Analysis: "
- << setw(5)
- << 100. * static_cast<double>(csnp)
- / static_cast<double>(nsnps)
- << "%...";
- } else
+ << setw(5)
+ << 100. * static_cast<double>(csnp)
+ / static_cast<double>(nsnps)
+ << "%...";
+ }
+ else
{
-
std::cout << "\b\b\b\b\b\b\b\b\b"
- << setw(5)
- << 100. * static_cast<double>(csnp)
- / static_cast<double>(nsnps)
- << "%...";
+ << setw(5)
+ << 100. * static_cast<double>(csnp)
+ / static_cast<double>(nsnps)
+ << "%...";
}
std::cout.flush();
}
@@ -78,7 +78,7 @@
}
void open_files_for_output(std::vector<std::ofstream*>& outfile,
- std::string& outfilename_str)
+ std::string& outfilename_str)
{
//create a list of filenames
const int amount_of_files = 5;
@@ -96,8 +96,8 @@
if (!outfile[i]->is_open())
{
std::cerr << "Cannot open file for writing: "
- << filenames[i]
- << "\n";
+ << filenames[i]
+ << "\n";
exit(1);
}
}
@@ -106,22 +106,25 @@
int create_phenotype(phedata& phd, cmdvars& input_var)
{
phd.set_is_interaction_excluded(input_var.isIsInteractionExcluded());
- phd.setphedata(input_var.getPhefilename(), input_var.getNoutcomes(),
- input_var.getNpeople(), input_var.getInteraction(),
- input_var.isIscox());
+ phd.setphedata(input_var.getPhefilename(),
+ input_var.getNoutcomes(),
+ input_var.getNpeople(),
+ input_var.getInteraction(),
+ input_var.isIscox());
int interaction_cox = input_var.getInteraction();
#if COXPH
interaction_cox--;
#endif
- if (input_var.getInteraction() < 0 || input_var.getInteraction() > phd.ncov
- || interaction_cox > phd.ncov)
+ if (input_var.getInteraction() < 0 ||
+ input_var.getInteraction() > phd.ncov ||
+ interaction_cox > phd.ncov)
{
std::cerr << "error: Interaction parameter is out of range "
- << "(interaction="
- << input_var.getInteraction()
- << ") \n";
+ << "(interaction="
+ << input_var.getInteraction()
+ << ") \n";
exit(1);
}
@@ -144,22 +147,22 @@
for (unsigned int i = 0; i < outfile.size(); i++)
{
(*outfile[i]) << "name"
- << input_var.getSep()
- << "A1"
- << input_var.getSep()
- << "A2"
- << input_var.getSep()
- << "Freq1"
- << input_var.getSep()
- << "MAF"
- << input_var.getSep()
- << "Quality"
- << input_var.getSep()
- << "Rsq"
- << input_var.getSep()
- << "n"
- << input_var.getSep()
- << "Mean_predictor_allele";
+ << input_var.getSep()
+ << "A1"
+ << input_var.getSep()
+ << "A2"
+ << input_var.getSep()
+ << "Freq1"
+ << input_var.getSep()
+ << "MAF"
+ << input_var.getSep()
+ << "Quality"
+ << input_var.getSep()
+ << "Rsq"
+ << input_var.getSep()
+ << "n"
+ << input_var.getSep()
+ << "Mean_predictor_allele";
if (input_var.getChrom() != "-1")
(*outfile[i]) << input_var.getSep() << "chrom";
if (input_var.getMapfilename() != NULL)
@@ -171,11 +174,11 @@
for (unsigned int file = 0; file < outfile.size(); file++)
for (int i = 0; i < phd.n_model_terms - 1; i++)
*outfile[file] << input_var.getSep()
- << "beta_"
- << phd.model_terms[i]
- << input_var.getSep()
- << "sebeta_"
- << phd.model_terms[i];
+ << "beta_"
+ << phd.model_terms[i]
+ << input_var.getSep()
+ << "sebeta_"
+ << phd.model_terms[i];
}
}
@@ -185,73 +188,72 @@
create_start_of_header(outfile, input_var, phd);
*outfile[0] << input_var.getSep()
- << "beta_SNP_A1A2"
- << input_var.getSep()
- << "beta_SNP_A1A1"
- << input_var.getSep()
- << "sebeta_SNP_A1A2"
- << input_var.getSep()
- << "sebeta_SNP_A1A1";
-
+ << "beta_SNP_A1A2"
+ << input_var.getSep()
+ << "beta_SNP_A1A1"
+ << input_var.getSep()
+ << "sebeta_SNP_A1A2"
+ << input_var.getSep()
+ << "sebeta_SNP_A1A1";
*outfile[1] << input_var.getSep()
- << "beta_SNP_addA1"
- << input_var.getSep()
- << "sebeta_SNP_addA1";
+ << "beta_SNP_addA1"
+ << input_var.getSep()
+ << "sebeta_SNP_addA1";
*outfile[2] << input_var.getSep()
- << "beta_SNP_domA1"
- << input_var.getSep()
- << "sebeta_SNP_domA1";
+ << "beta_SNP_domA1"
+ << input_var.getSep()
+ << "sebeta_SNP_domA1";
*outfile[3] << input_var.getSep()
- << "beta_SNP_recA1"
- << input_var.getSep()
- << "sebeta_SNP_recA1";
+ << "beta_SNP_recA1"
+ << input_var.getSep()
+ << "sebeta_SNP_recA1";
*outfile[4] << input_var.getSep()
- << "beta_SNP_odom"
- << input_var.getSep()
- << "sebeta_SNP_odom";
+ << "beta_SNP_odom"
+ << input_var.getSep()
+ << "sebeta_SNP_odom";
if (input_var.getInteraction() != 0)
{
//Han Chen
*outfile[0] << input_var.getSep()
- << "beta_SNP_A1A2_"
- << phd.model_terms[interaction_cox]
- << input_var.getSep()
- << "sebeta_SNP_A1A2_"
- << phd.model_terms[interaction_cox]
- << input_var.getSep()
- << "beta_SNP_A1A1_"
- << phd.model_terms[interaction_cox]
- << input_var.getSep()
- << "sebeta_SNP_A1A1_"
- << phd.model_terms[interaction_cox];
+ << "beta_SNP_A1A2_"
+ << phd.model_terms[interaction_cox]
+ << input_var.getSep()
+ << "sebeta_SNP_A1A2_"
+ << phd.model_terms[interaction_cox]
+ << input_var.getSep()
+ << "beta_SNP_A1A1_"
+ << phd.model_terms[interaction_cox]
+ << input_var.getSep()
+ << "sebeta_SNP_A1A1_"
+ << phd.model_terms[interaction_cox];
#if !COXPH
if (input_var.getInverseFilename() == NULL && !input_var.getAllcov())
{
*outfile[0] << input_var.getSep()
- << "cov_SNP_A1A2_int_SNP_"
- << phd.model_terms[interaction_cox]
- << input_var.getSep()
- << "cov_SNP_A1A1_int_SNP_"
- << phd.model_terms[interaction_cox];
+ << "cov_SNP_A1A2_int_SNP_"
+ << phd.model_terms[interaction_cox]
+ << input_var.getSep()
+ << "cov_SNP_A1A1_int_SNP_"
+ << phd.model_terms[interaction_cox];
}
#endif
//Oct 26, 2009
for (unsigned int file = 1; file < outfile.size(); file++)
{
*outfile[file] << input_var.getSep()
- << "beta_SNP_"
- << phd.model_terms[interaction_cox]
- << input_var.getSep()
- << "sebeta_SNP_"
- << phd.model_terms[interaction_cox];
+ << "beta_SNP_"
+ << phd.model_terms[interaction_cox]
+ << input_var.getSep()
+ << "sebeta_SNP_"
+ << phd.model_terms[interaction_cox];
//Han Chen
#if !COXPH
if (input_var.getInverseFilename() == NULL
&& !input_var.getAllcov())
{
*outfile[file] << input_var.getSep()
- << "cov_SNP_int_SNP_"
- << phd.model_terms[interaction_cox];
+ << "cov_SNP_int_SNP_"
+ << phd.model_terms[interaction_cox];
}
#endif
//Oct 26, 2009
@@ -265,22 +267,22 @@
}
void create_header2(std::vector<std::ofstream*>& outfile, cmdvars& input_var,
- phedata& phd, int interaction_cox)
+ phedata& phd, int interaction_cox)
{
create_start_of_header(outfile, input_var, phd);
*outfile[0] << input_var.getSep()
- << "beta_SNP_add"
- << input_var.getSep()
- << "sebeta_SNP_add";
+ << "beta_SNP_add"
+ << input_var.getSep()
+ << "sebeta_SNP_add";
if (input_var.getInteraction() != 0)
{
*outfile[0] << input_var.getSep()
- << "beta_SNP_"
- << phd.model_terms[interaction_cox]
- << input_var.getSep()
- << "sebeta_SNP_"
- << phd.model_terms[interaction_cox];
+ << "beta_SNP_"
+ << phd.model_terms[interaction_cox]
+ << input_var.getSep()
+ << "sebeta_SNP_"
+ << phd.model_terms[interaction_cox];
}
if (input_var.getInverseFilename() == NULL)
@@ -290,8 +292,8 @@
if (input_var.getInteraction() != 0 && !input_var.getAllcov())
{
*outfile[0] << input_var.getSep()
- << "cov_SNP_int_SNP_"
- << phd.model_terms[interaction_cox];
+ << "cov_SNP_int_SNP_"
+ << phd.model_terms[interaction_cox];
}
#endif
*outfile[0] << input_var.getSep() << "loglik"; //"chi2_SNP";
@@ -301,26 +303,26 @@
}
void write_mlinfo(const std::vector<std::ofstream*>& outfile, unsigned int file,
- const mlinfo& mli, int csnp, const cmdvars& input_var, int gcount,
- double freq)
+ const mlinfo& mli, int csnp, const cmdvars& input_var,
+ int gcount, double freq)
{
*outfile[file] << mli.name[csnp]
- << input_var.getSep()
- << mli.A1[csnp]
- << input_var.getSep()
- << mli.A2[csnp]
- << input_var.getSep()
- << mli.Freq1[csnp]
- << input_var.getSep()
- << mli.MAF[csnp]
- << input_var.getSep()
- << mli.Quality[csnp]
- << input_var.getSep()
- << mli.Rsq[csnp]
- << input_var.getSep()
- << gcount
- << input_var.getSep()
- << freq;
+ << input_var.getSep()
+ << mli.A1[csnp]
+ << input_var.getSep()
+ << mli.A2[csnp]
+ << input_var.getSep()
+ << mli.Freq1[csnp]
+ << input_var.getSep()
+ << mli.MAF[csnp]
+ << input_var.getSep()
+ << mli.Quality[csnp]
+ << input_var.getSep()
+ << mli.Rsq[csnp]
+ << input_var.getSep()
+ << gcount
+ << input_var.getSep()
+ << freq;
if (input_var.getChrom() != "-1")
{
*outfile[file] << input_var.getSep() << input_var.getChrom();
@@ -335,7 +337,9 @@
int number_of_rows_or_columns)
{
int start_pos;
- if (!input_var.getAllcov() && model == 0 && input_var.getInteraction() == 0)
+ if (!input_var.getAllcov() &&
+ model == 0 &&
+ input_var.getInteraction() == 0)
{
if (input_var.getNgpreds() == 2)
{
@@ -382,11 +386,8 @@
phedata phd;
int interaction_cox = create_phenotype(phd, input_var);
-
masked_matrix invvarmatrix;
-
-
std::cout << "Reading data ..." << std::flush;
if (input_var.getInverseFilename() != NULL)
{
@@ -412,7 +413,6 @@
std::cout << " loaded genotypic data ..." << std::flush;
-
// estimate null model
#if COXPH
coxph_data nrgd = coxph_data(phd, gtd, -1);
@@ -436,12 +436,12 @@
std::cout << "[DEBUG] linear_reg nrd = linear_reg(nrgd); DONE.";
#endif
nrd.estimate(nrgd, 0, CHOLTOL, 0, input_var.getInteraction(),
- input_var.getNgpreds(), invvarmatrix, input_var.getRobust(), 1);
+ input_var.getNgpreds(), invvarmatrix,
+ input_var.getRobust(), 1);
#elif COXPH
- coxph_reg nrd(nrgd);
-
+ coxph_reg nrd = coxph_reg(nrgd);
nrd.estimate(nrgd, 0, MAXITER, EPS, CHOLTOL, 0,
- input_var.getInteraction(), input_var.getNgpreds(), 1);
+ input_var.getInteraction(), input_var.getNgpreds(), 1);
#endif
std::cout << " estimated null model ...";
@@ -451,7 +451,6 @@
#else
regdata rgd(phd, gtd, 0, input_var.isIsInteractionExcluded());
#endif
-
std::cout << " formed regression object ...";
//________________________________________________________________
@@ -470,13 +469,13 @@
} else //Only additive model. Only one output file
{
outfile.push_back(
- new std::ofstream((outfilename_str + "_add.out.txt").c_str()));
+ new std::ofstream((outfilename_str + "_add.out.txt").c_str()));
if (!outfile[0]->is_open())
{
std::cerr << "Cannot open file for writing: "
- << outfilename_str
- << "\n";
+ << outfilename_str
+ << "\n";
exit(1);
}
if (input_var.getNohead() != 1)
@@ -534,11 +533,13 @@
// freq = (gtd.G).column_mean(csnp)/2.;
gtd.get_var(csnp, snpdata1);
for (unsigned int ii = 0; ii < gtd.nids; ii++)
+ {
if (!isnan(snpdata1[ii]))
{
gcount++;
freq += snpdata1[ii] * 0.5;
}
+ }
}
freq /= static_cast<double>(gcount);
int poly = 1;
@@ -559,7 +560,8 @@
// Write mlinfo to output:
for (unsigned int file = 0; file < outfile.size(); file++)
{
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/genabel -r 1231
More information about the Genabel-commits
mailing list