[Genabel-commits] r1230 - in branches/ProbABEL-pacox/v.0.3.0/ProbABEL: . examples src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sun May 19 23:27:46 CEST 2013
Author: lckarssen
Date: 2013-05-19 23:27:46 +0200 (Sun, 19 May 2013)
New Revision: 1230
Modified:
branches/ProbABEL-pacox/v.0.3.0/ProbABEL/configure.ac
branches/ProbABEL-pacox/v.0.3.0/ProbABEL/examples/example_mms.sh
branches/ProbABEL-pacox/v.0.3.0/ProbABEL/src/Makefile.am
branches/ProbABEL-pacox/v.0.3.0/ProbABEL/src/command_line_settings.cpp
branches/ProbABEL-pacox/v.0.3.0/ProbABEL/src/command_line_settings.h
branches/ProbABEL-pacox/v.0.3.0/ProbABEL/src/coxph_data.cpp
branches/ProbABEL-pacox/v.0.3.0/ProbABEL/src/eigen_mematrix.cpp
branches/ProbABEL-pacox/v.0.3.0/ProbABEL/src/gendata.cpp
branches/ProbABEL-pacox/v.0.3.0/ProbABEL/src/gendata.h
branches/ProbABEL-pacox/v.0.3.0/ProbABEL/src/main.cpp
branches/ProbABEL-pacox/v.0.3.0/ProbABEL/src/mematri1.h
branches/ProbABEL-pacox/v.0.3.0/ProbABEL/src/prepare_data.R
branches/ProbABEL-pacox/v.0.3.0/ProbABEL/src/reg1.cpp
branches/ProbABEL-pacox/v.0.3.0/ProbABEL/src/reg1.h
branches/ProbABEL-pacox/v.0.3.0/ProbABEL/src/regdata.cpp
Log:
ProbABEL pacoxfix v0.3.0 branch: merge changes from trunk.
Modified: branches/ProbABEL-pacox/v.0.3.0/ProbABEL/configure.ac
===================================================================
--- branches/ProbABEL-pacox/v.0.3.0/ProbABEL/configure.ac 2013-05-19 20:00:36 UTC (rev 1229)
+++ branches/ProbABEL-pacox/v.0.3.0/ProbABEL/configure.ac 2013-05-19 21:27:46 UTC (rev 1230)
@@ -2,7 +2,7 @@
# Process this file with autoconf to produce a configure script.
AC_PREREQ([2.67])
-AC_INIT(ProbABEL, 0.3.0-pacox-beta-1, genabel-devel at r-forge.wu-wien.ac.at)
+AC_INIT(ProbABEL, 0.3.0, genabel-devel at r-forge.wu-wien.ac.at)
AM_INIT_AUTOMAKE([silent-rules])
AM_SILENT_RULES
AC_CONFIG_SRCDIR([src/data.h])
@@ -40,29 +40,28 @@
# Checks for header files.
AC_FUNC_ALLOCA
AC_CHECK_HEADERS([float.h inttypes.h libintl.h limits.h stddef.h \
- stdint.h stdlib.h string.h sys/param.h wchar.h wctype.h])
+ stdint.h stdlib.h string.h sys/param.h wchar.h wctype.h])
# See if we want use of the Eigen library enabled (yes by default) and if so,
# whether we can find the library files.
AC_ARG_WITH([eigen],
- [AS_HELP_STRING([--with-eigen], [Use the Eigen template library for fast \
- linear algebra (Eigen can be downloaded \
- from eigen.tuxfamily.org); this is enabled \
- by default])],
- [eigen=${enableval}],
- [eigen=yes])
+ AS_HELP_STRING([--with-eigen], [Use the Eigen template library for fast \
+ linear algebra (Eigen can be downloaded \
+ from eigen.tuxfamily.org); this is enabled \
+ by default]
+ )
+)
-if test "x$eigen" = "xyes"; then
+if test "x$with_eigen" != "xno"; then
AC_MSG_NOTICE([building using the Eigen headers enabled])
AC_ARG_WITH([eigen-include-path],
- [AS_HELP_STRING([--with-eigen-include-path],
- [location of the Eigen headers, defaults to /usr/include/eigen3])],
- [CXXFLAGS+=" -I${withval}"],
- [CXXFLAGS+=' -I/usr/include/eigen3'
- CPPFLAGS+=' -I/usr/include/eigen3'])
- dnl AC_SUBST([EIGEN_CFLAGS])
+ [AS_HELP_STRING([--with-eigen-include-path],
+ [location of the Eigen headers, defaults to /usr/include/eigen3])],
+ [CXXFLAGS+=" -I${withval}"],
+ [CXXFLAGS+=' -I/usr/include/eigen3'
+ CPPFLAGS+=' -I/usr/include/eigen3'])
# Check for the EIGEN header files
AC_CHECK_HEADERS([Eigen/Dense Eigen/LU])
@@ -75,7 +74,7 @@
else
AC_MSG_NOTICE([not using Eigen for linear algebra])
fi
-AM_CONDITIONAL([WITH_EIGEN], test "x$eigen" = "xyes")
+AM_CONDITIONAL([WITH_EIGEN], test "x$with_eigen" != "xno")
# Checks for typedefs, structures, and compiler characteristics.
@@ -124,11 +123,11 @@
# Files to be generated by autotools
AC_CONFIG_FILES([
- Makefile
- src/Makefile
- doc/Makefile
- examples/Makefile
- tests/Makefile
+ Makefile
+ src/Makefile
+ doc/Makefile
+ examples/Makefile
+ tests/Makefile
])
# Create output files
AC_OUTPUT
Modified: branches/ProbABEL-pacox/v.0.3.0/ProbABEL/examples/example_mms.sh
===================================================================
--- branches/ProbABEL-pacox/v.0.3.0/ProbABEL/examples/example_mms.sh 2013-05-19 20:00:36 UTC (rev 1229)
+++ branches/ProbABEL-pacox/v.0.3.0/ProbABEL/examples/example_mms.sh 2013-05-19 21:27:46 UTC (rev 1230)
@@ -35,8 +35,8 @@
>& 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"
@@ -60,8 +60,8 @@
for model in add domin over_domin recess 2df; do
run_diff mmscore_prob_${model}.out.txt \
- mmscore_prob_fv_${model}.out.txt \
- "mmscore check ($model model): prob vs. prob_fv"
+ mmscore_prob_fv_${model}.out.txt \
+ "mmscore check ($model model): prob vs. prob_fv"
done
# The following checks are disabled because of the missing LogLik
Modified: branches/ProbABEL-pacox/v.0.3.0/ProbABEL/src/Makefile.am
===================================================================
--- branches/ProbABEL-pacox/v.0.3.0/ProbABEL/src/Makefile.am 2013-05-19 20:00:36 UTC (rev 1229)
+++ branches/ProbABEL-pacox/v.0.3.0/ProbABEL/src/Makefile.am 2013-05-19 21:27:46 UTC (rev 1230)
@@ -6,7 +6,7 @@
REGFILES = data.h data.cpp gendata.h gendata.cpp mematrix.h mematri1.h \
command_line_settings.h command_line_settings.cpp reg1.h usage.h \
usage.cpp main.cpp utilities.h utilities.cpp phedata.h phedata.cpp \
- cholesky.h cholesky.cpp regdata.h regdata.cpp \
+ cholesky.h cholesky.cpp regdata.h regdata.cpp \
maskedmatrix.cpp maskedmatrix.h reg1.cpp
EIGENFILES = eigen_mematrix.h eigen_mematrix.cpp
@@ -70,6 +70,10 @@
## Install these scripts in the bin directory as well:
dist_bin_SCRIPTS = probabel.pl extIDS.pl
+## Install this R script in the examples directory
+scriptdir = $(pkgdatadir)/scripts
+dist_script_DATA = prepare_data.R
+
## Install the config file
dist_sysconf_DATA = probabel_config.cfg.example
Modified: branches/ProbABEL-pacox/v.0.3.0/ProbABEL/src/command_line_settings.cpp
===================================================================
--- branches/ProbABEL-pacox/v.0.3.0/ProbABEL/src/command_line_settings.cpp 2013-05-19 20:00:36 UTC (rev 1229)
+++ branches/ProbABEL-pacox/v.0.3.0/ProbABEL/src/command_line_settings.cpp 2013-05-19 21:27:46 UTC (rev 1230)
@@ -93,7 +93,7 @@
return npeople;
}
-char* cmdvars::getOutfilename() const
+string cmdvars::getOutfilename() const
{
return outfilename;
}
@@ -227,7 +227,7 @@
break;
case '?':
- print_usage(program_name, 1);
+ print_usage(program_name, 2);
case -1:
break;
default:
@@ -248,13 +248,30 @@
<< "\n(C) Yurii Aulchenko, Lennart C. Karssen, Maksim Struchalin, "
<< "EMCR\n\n";
#if EIGEN
- cout << "Using EIGEN version "<<EIGEN_WORLD_VERSION
- <<"."<<EIGEN_MAJOR_VERSION<<"."<<EIGEN_MINOR_VERSION
+ cout << "Using EIGEN version " << EIGEN_WORLD_VERSION
+ << "." << EIGEN_MAJOR_VERSION << "." << EIGEN_MINOR_VERSION
<< " for matrix operations\n";
#endif
if (neco[0] != 1 || neco[1] != 1 || neco[2] != 1)
{
+ cerr << endl;
+ if (neco[0] != 1)
+ {
+ cerr << "Error: Missing required phenotype file (-p/--pheno option)"
+ << endl;
+ }
+ if (neco[1] != 1)
+ {
+ cerr << "Error: Missing required info file (-i/--info option)"
+ << endl;
+ }
+ if (neco[2] != 1)
+ {
+ cerr << "Error: Missing required genotype file (-d/--dose option)"
+ << endl;
+ }
+ cerr << endl;
print_usage(program_name, 1);
}
@@ -298,7 +315,7 @@
cout << "\t --chrom = " << chrom << endl;
else
cout << "\t --chrom = not in output\n";
- if (outfilename != NULL)
+ if (outfilename.compare("") != 0)
cout << "\t --out = " << outfilename << endl;
else
cout << "\t --out = " << "regression.out.txt" << endl;
@@ -333,9 +350,9 @@
interaction = interaction_excluded; //ups
is_interaction_excluded = true;
}
- if (outfilename == NULL)
+ if (outfilename.compare("") == 0)
{
- outfilename = (char *) string("regression").c_str();
+ outfilename = string("regression");
}
#if COXPH
if (score)
Modified: branches/ProbABEL-pacox/v.0.3.0/ProbABEL/src/command_line_settings.h
===================================================================
--- branches/ProbABEL-pacox/v.0.3.0/ProbABEL/src/command_line_settings.h 2013-05-19 20:00:36 UTC (rev 1229)
+++ branches/ProbABEL-pacox/v.0.3.0/ProbABEL/src/command_line_settings.h 2013-05-19 21:27:46 UTC (rev 1230)
@@ -20,7 +20,7 @@
char *mlinfofilename;
char *genfilename;
char *mapfilename;
- char *outfilename;
+ string outfilename;
char *inverse_filename;
string str_genfilename;
@@ -35,7 +35,7 @@
int robust;
string chrom;
string sep;
- int neco[3];
+ int neco[3]; /* Necessary command line options */
bool iscox;
int isFVF;
int noutcomes;
@@ -52,7 +52,7 @@
mlinfofilename = NULL;
genfilename = NULL;
mapfilename = NULL;
- outfilename = NULL;
+ outfilename = string("");
inverse_filename = NULL;
sep = " ";
@@ -94,7 +94,7 @@
int getNohead() const;
int getNoutcomes() const;
int getNpeople() const;
- char* getOutfilename() const;
+ string getOutfilename() const;
char* getProgramName() const;
int getRobust() const;
int getScore() const;
Modified: branches/ProbABEL-pacox/v.0.3.0/ProbABEL/src/coxph_data.cpp
===================================================================
--- branches/ProbABEL-pacox/v.0.3.0/ProbABEL/src/coxph_data.cpp 2013-05-19 20:00:36 UTC (rev 1229)
+++ branches/ProbABEL-pacox/v.0.3.0/ProbABEL/src/coxph_data.cpp 2013-05-19 21:27:46 UTC (rev 1230)
@@ -114,7 +114,7 @@
{
for (int j = 0; j < ngpreds; j++)
{
- float snpdata[nids];
+ double snpdata[nids];
gend.get_var(snpnum * ngpreds + j, snpdata);
for (int i = 0; i < nids; i++)
X.put(snpdata[i], i, (ncov - ngpreds + j));
@@ -195,7 +195,7 @@
for (int j = 0; j < ngpreds; j++)
{
- float snpdata[nids];
+ double snpdata[nids];
for (int i = 0; i < nids; i++)
{
masked_data[i] = 0;
Modified: branches/ProbABEL-pacox/v.0.3.0/ProbABEL/src/eigen_mematrix.cpp
===================================================================
--- branches/ProbABEL-pacox/v.0.3.0/ProbABEL/src/eigen_mematrix.cpp 2013-05-19 20:00:36 UTC (rev 1229)
+++ branches/ProbABEL-pacox/v.0.3.0/ProbABEL/src/eigen_mematrix.cpp 2013-05-19 21:27:46 UTC (rev 1230)
@@ -273,8 +273,9 @@
{
std:: cout << "nr=" << i << ":\t";
for (int j = 0; j < ncol; j++)
-// cout << data.data()[i * ncol + j] << "\t";
- printf("%f\t", data.data()[i * ncol + j]);
+ {
+ printf("%1.9e\t", data.data()[i * ncol + j]);
+ }
std::cout << "\n";
}
}
Modified: branches/ProbABEL-pacox/v.0.3.0/ProbABEL/src/gendata.cpp
===================================================================
--- branches/ProbABEL-pacox/v.0.3.0/ProbABEL/src/gendata.cpp 2013-05-19 20:00:36 UTC (rev 1229)
+++ branches/ProbABEL-pacox/v.0.3.0/ProbABEL/src/gendata.cpp 2013-05-19 21:27:46 UTC (rev 1230)
@@ -5,6 +5,7 @@
* Author: mkooyman
*/
#include <string>
+#include <errno.h>
#include "gendata.h"
#include "fvlib/FileVector.h"
#if EIGEN
@@ -16,29 +17,53 @@
#endif
#include "utilities.h"
-void gendata::get_var(int var, float * data)
+void gendata::get_var(int var, double * data)
{
- if (DAG == NULL)
+ if (DAG == NULL) // Read from text file
{
for (int i = 0; i < G.nrow; i++)
{
data[i] = G.get(i, var);
}
}
- else if (DAG != NULL)
+ else if (DAG != NULL) // Read from fv file
{
- float tmpdata[DAG->getNumObservations()];
+ double tmpdata[DAG->getNumObservations()];
DAG->readVariableAs((unsigned long int) var, tmpdata);
+
unsigned int j = 0;
for (unsigned int i = 0; i < DAG->getNumObservations(); i++)
{
if (!DAGmask[i])
{
- data[j++] = tmpdata[i];
+ // A dirty trick to get rid of conversion
+ // errors. Instead of casting float data to double we
+ // convert the data to string and then do strtod()
+ std::ostringstream strs;
+ strs << tmpdata[i];
+ std::string str = strs.str();
+ double val;
+ char *endptr;
+ errno = 0; // To distinguish success/failure
+ // after strtod()
+ val = strtod(str.c_str(), &endptr);
+
+ if ((errno == ERANGE && (val == HUGE_VALF || val == HUGE_VALL))
+ || (errno != 0 && val == 0)) {
+ perror("Error while reading genetic data (strtod)");
+ exit(EXIT_FAILURE);
+ }
+
+ if (endptr == str.c_str()) {
+ cerr << "No digits were found while reading genetic data"
+ << endl;
+ exit(EXIT_FAILURE);
+ }
+
+ /* If we got here, strtod() successfully parsed a number */
+ data[j++] = val;
}
}
- // std::cout << j << " " << DAG->getNumObservations() << " "
- // << nids << "\n";
}
else
{
@@ -165,10 +190,27 @@
{
infile >> tmpstr;
// tmpstr contains the dosage in string form. Convert
- // it to float (if tmpstr is NaN it will be set to nan).
- // Note that Valgrind 3.7.0 gives "Invalid read of
- // size 8" error messages here. A bug in Valgrind?!
- float dosage = strtod(tmpstr.c_str(), (char **) NULL);
+ // 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 ((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"
+ << endl;
+ exit(EXIT_FAILURE);
+ }
+
+ /* If we got here, strtod() successfully parsed a number */
G.put(dosage, k, j);
}
else
Modified: branches/ProbABEL-pacox/v.0.3.0/ProbABEL/src/gendata.h
===================================================================
--- branches/ProbABEL-pacox/v.0.3.0/ProbABEL/src/gendata.h 2013-05-19 20:00:36 UTC (rev 1229)
+++ branches/ProbABEL-pacox/v.0.3.0/ProbABEL/src/gendata.h 2013-05-19 21:27:46 UTC (rev 1230)
@@ -32,7 +32,7 @@
unsigned int npeople, unsigned int nmeasured,
unsigned short int * allmeasured, std::string * idnames);
- void get_var(int var, float * data);
+ void get_var(int var, double * data);
~gendata();
@@ -40,7 +40,7 @@
// ANOTHER PRIVATE OBJECT IS A POINTER TO DATABELBASECPP
// UPDATE SNP, ALL REGRESSION METHODS: ACCOUNT FOR MISSING
private:
- mematrix<float> G;
+ mematrix<double> G;
AbstractMatrix * DAG;
unsigned short int * DAGmask;
};
Modified: branches/ProbABEL-pacox/v.0.3.0/ProbABEL/src/main.cpp
===================================================================
--- branches/ProbABEL-pacox/v.0.3.0/ProbABEL/src/main.cpp 2013-05-19 20:00:36 UTC (rev 1229)
+++ branches/ProbABEL-pacox/v.0.3.0/ProbABEL/src/main.cpp 2013-05-19 21:27:46 UTC (rev 1230)
@@ -59,19 +59,22 @@
if (csnp == 0)
{
std::cout << "Analysis: "
- << setw(5) << 100. * static_cast<double>(csnp)
+ << 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)
+ << setw(5)
+ << 100. * static_cast<double>(csnp)
/ static_cast<double>(nsnps)
<< "%...";
}
std::cout.flush();
}
+ std::cout << setprecision(6);
}
void open_files_for_output(std::vector<std::ofstream*>& outfile,
@@ -79,12 +82,12 @@
{
//create a list of filenames
const int amount_of_files = 5;
- std::string filenames[amount_of_files] =
- { outfilename_str + "_2df.out.txt",
- outfilename_str + "_add.out.txt",
- outfilename_str + "_domin.out.txt",
- outfilename_str + "_recess.out.txt",
- outfilename_str + "_over_domin.out.txt" };
+ std::string filenames[amount_of_files] = {
+ outfilename_str + "_2df.out.txt",
+ outfilename_str + "_add.out.txt",
+ outfilename_str + "_domin.out.txt",
+ outfilename_str + "_recess.out.txt",
+ outfilename_str + "_over_domin.out.txt" };
for (int i = 0; i < amount_of_files; i++)
{
@@ -92,7 +95,8 @@
outfile[i]->open((filenames[i]).c_str());
if (!outfile[i]->is_open())
{
- std::cerr << "Cannot open file for writing: " << filenames[i]
+ std::cerr << "Cannot open file for writing: "
+ << filenames[i]
<< "\n";
exit(1);
}
@@ -118,7 +122,9 @@
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);
}
@@ -140,14 +146,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()
+ (*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";
if (input_var.getChrom() != "-1")
(*outfile[i]) << input_var.getSep() << "chrom";
@@ -160,8 +174,10 @@
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_"
+ << "beta_"
+ << phd.model_terms[i]
+ << input_var.getSep()
+ << "sebeta_"
<< phd.model_terms[i];
}
}
@@ -171,51 +187,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";
-
- *outfile[1] << input_var.getSep() << "beta_SNP_addA1" << input_var.getSep()
+ *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";
+ *outfile[1] << input_var.getSep()
+ << "beta_SNP_addA1"
+ << input_var.getSep()
<< "sebeta_SNP_addA1";
- *outfile[2] << input_var.getSep() << "beta_SNP_domA1" << input_var.getSep()
+ *outfile[2] << input_var.getSep()
+ << "beta_SNP_domA1"
+ << input_var.getSep()
<< "sebeta_SNP_domA1";
- *outfile[3] << input_var.getSep() << "beta_SNP_recA1" << input_var.getSep()
+ *outfile[3] << input_var.getSep()
+ << "beta_SNP_recA1"
+ << input_var.getSep()
<< "sebeta_SNP_recA1";
- *outfile[4] << input_var.getSep() << "beta_SNP_odom" << input_var.getSep()
+ *outfile[4] << input_var.getSep()
+ << "beta_SNP_odom"
+ << input_var.getSep()
<< "sebeta_SNP_odom";
if (input_var.getInteraction() != 0)
{
//Han Chen
- *outfile[0] << input_var.getSep() << "beta_SNP_A1A2_"
+ *outfile[0] << input_var.getSep()
+ << "beta_SNP_A1A2_"
<< phd.model_terms[interaction_cox]
- << input_var.getSep() << "sebeta_SNP_A1A2_"
+ << input_var.getSep()
+ << "sebeta_SNP_A1A2_"
<< phd.model_terms[interaction_cox]
- << input_var.getSep() << "beta_SNP_A1A1_"
+ << input_var.getSep()
+ << "beta_SNP_A1A1_"
<< phd.model_terms[interaction_cox]
- << input_var.getSep() << "sebeta_SNP_A1A1_"
+ << 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_"
+ *outfile[0] << input_var.getSep()
+ << "cov_SNP_A1A2_int_SNP_"
<< phd.model_terms[interaction_cox]
- << input_var.getSep() << "cov_SNP_A1A1_int_SNP_"
+ << 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_"
+ *outfile[file] << input_var.getSep()
+ << "beta_SNP_"
<< phd.model_terms[interaction_cox]
- << input_var.getSep() << "sebeta_SNP_"
+ << input_var.getSep()
+ << "sebeta_SNP_"
<< phd.model_terms[interaction_cox];
//Han Chen
#if !COXPH
if (input_var.getInverseFilename() == NULL
- && !input_var.getAllcov())
+ && !input_var.getAllcov())
{
- *outfile[file] << input_var.getSep() << "cov_SNP_int_SNP_"
+ *outfile[file] << input_var.getSep()
+ << "cov_SNP_int_SNP_"
<< phd.model_terms[interaction_cox];
}
#endif
@@ -233,14 +270,18 @@
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";
+ *outfile[0] << input_var.getSep()
+ << "beta_SNP_add"
+ << input_var.getSep()
+ << "sebeta_SNP_add";
if (input_var.getInteraction() != 0)
{
- *outfile[0] << input_var.getSep() << "beta_SNP_"
+ *outfile[0] << input_var.getSep()
+ << "beta_SNP_"
<< phd.model_terms[interaction_cox]
- << input_var.getSep() << "sebeta_SNP_"
+ << input_var.getSep()
+ << "sebeta_SNP_"
<< phd.model_terms[interaction_cox];
}
@@ -250,7 +291,8 @@
#if !COXPH
if (input_var.getInteraction() != 0 && !input_var.getAllcov())
{
- *outfile[0] << input_var.getSep() << "cov_SNP_int_SNP_"
+ *outfile[0] << input_var.getSep()
+ << "cov_SNP_int_SNP_"
<< phd.model_terms[interaction_cox];
}
#endif
@@ -260,53 +302,92 @@
*outfile[0] << "\n";
}
+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)
+{
+ *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;
+ if (input_var.getChrom() != "-1")
+ {
+ *outfile[file] << input_var.getSep() << input_var.getChrom();
+ }
+ if (input_var.getMapfilename() != NULL)
+ {
+ *outfile[file] << input_var.getSep() << mli.map[csnp];
+ }
+}
+
+int get_start_position(const cmdvars& input_var, int model,
+ int number_of_rows_or_columns)
+{
+ int start_pos;
+ if (!input_var.getAllcov() &&
+ model == 0 &&
+ input_var.getInteraction() == 0)
+ {
+ if (input_var.getNgpreds() == 2)
+ {
+ start_pos = number_of_rows_or_columns - 2;
+ } else
+ {
+ start_pos = number_of_rows_or_columns - 1;
+ }
+ } else if (!input_var.getAllcov() && model == 0
+ && input_var.getInteraction() != 0)
+ {
+ if (input_var.getNgpreds() == 2)
+ {
+ start_pos = number_of_rows_or_columns - 4;
+ } else
+ {
+ start_pos = number_of_rows_or_columns - 2;
+ }
+ } else if (!input_var.getAllcov() && model != 0
+ && input_var.getInteraction() == 0)
+ {
+ start_pos = number_of_rows_or_columns - 1;
+ } else if (!input_var.getAllcov() && model != 0
+ && input_var.getInteraction() != 0)
+ {
+ start_pos = number_of_rows_or_columns - 2;
+ } else
+ {
+ start_pos = 0;
+ }
+
+ return start_pos;
+}
+
int main(int argc, char * argv[])
{
cmdvars input_var;
input_var.set_variables(argc, argv);
input_var.printinfo();
- // if (allcov && ngpreds>1)
- // {
- // cout << "\n\n"
- // << "WARNING: --allcov allowed only for 1 predictor (MLDOSE)\n";
- // allcov = 0;
- // }
+
mlinfo mli(input_var.getMlinfofilename(), input_var.getMapfilename());
int nsnps = mli.nsnps;
phedata phd;
int interaction_cox = create_phenotype(phd, input_var);
- //interaction--;
- // if (input_var.getInverseFilename()!= NULL && phd.ncov > 1)
- // {
- // std::cerr << "Error: In mmscore you can not use any covariates."
- // << " You phenotype file must conatin id column and "
- // << "trait (residuals) only\n";
- // exit(1);
- // }
- // if (input_var.getInverseFilename()!= NULL &&
- // (allcov == 1 || score == 1
- // || input_var.getInteraction()!= 0
- // || ngpreds==2))
- // {
- // std::cerr << "Error: In mmscore you can use additive model "
- // << "without any inetractions only\n";
- // exit(1);
- // }
masked_matrix invvarmatrix;
- /*
- * now should be possible... delete this part later when everything works
- #if LOGISTIC
- if (input_var.getInverseFilename()!= NULL)
- {
- std::cerr << "ERROR: mmscore is forbidden for logistic regression\n";
- exit(1);
- }
- #endif
- */
-
std::cout << "Reading data ..." << std::flush;
if (input_var.getInverseFilename() != NULL)
{
@@ -315,27 +396,22 @@
gendata gtd;
if (!input_var.getIsFvf())
+ {
// use the non-filevector input format
gtd.re_gendata(input_var.getGenfilename(), nsnps,
input_var.getNgpreds(), phd.nids_all, phd.nids,
phd.allmeasured, input_var.getSkipd(), phd.idnames);
+ }
else
+ {
// use the filevector input format (missing second last skipd
// parameter)
gtd.re_gendata(input_var.getStrGenfilename(), nsnps,
input_var.getNgpreds(), phd.nids_all, phd.nids,
phd.allmeasured, phd.idnames);
+ }
std::cout << " loaded genotypic data ..." << std::flush;
- /**
- if (input_var.getIsFvf())
- gendata gtd(str_genfilename, nsnps, input_var.getNgpreds(),
- phd.nids_all, phd.allmeasured, phd.idnames);
- else
- gendata gtd(input_var.getGenfilename(), nsnps,
- input_var.getNgpreds(), phd.nids_all, phd.nids,
- phd.allmeasured, skipd, phd.idnames);
- **/
// estimate null model
#if COXPH
@@ -348,8 +424,11 @@
#if LOGISTIC
logistic_reg nrd = logistic_reg(nrgd);
nrd.estimate(nrgd, 0, MAXITER, EPS, CHOLTOL, 0,
- input_var.getInteraction(), input_var.getNgpreds(),
- invvarmatrix, input_var.getRobust(), 1);
+ input_var.getInteraction(),
+ input_var.getNgpreds(),
+ invvarmatrix,
+ input_var.getRobust(),
+ 1);
#elif LINEAR
linear_reg nrd = linear_reg(nrgd);
@@ -357,7 +436,8 @@
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 = coxph_reg(nrgd);
// LCK std::cout << "Starting estimation of null model...\n";
@@ -373,8 +453,6 @@
regdata rgd(phd, gtd, 0, input_var.isIsInteractionExcluded());
#endif
std::cout << " formed regression object ...";
-// rgd.X.print();
- std::cout << " done\n" << std::flush;
//________________________________________________________________
//Maksim, 9 Jan, 2009
@@ -389,15 +467,15 @@
{
create_header_1(outfile, input_var, phd, interaction_cox);
}
- }
- else //Only additive model. Only one output file
+ } else //Only additive model. Only one output file
{
outfile.push_back(
new std::ofstream((outfilename_str + "_add.out.txt").c_str()));
if (!outfile[0]->is_open())
{
- std::cerr << "Cannot open file for writing: " << outfilename_str
+ std::cerr << "Cannot open file for writing: "
+ << outfilename_str
<< "\n";
exit(1);
}
@@ -407,54 +485,7 @@
}
}
- //________________________________________________________________
- /*
- if (input_var.getAllcov())
- {
- if (score)
- {
- outfile << input_var.getSep() << "beta_mu"; // << input_var.getSep() << "beta_SNP_A1";
- outfile << input_var.getSep() << "sebeta_mu"; // << input_var.getSep() << "sebeta_SNP_A1";
- }
- else
- {
- for (int i =0; i<phd.n_model_terms-1;i++)
- outfile << input_var.getSep() << "beta_" << phd.model_terms[i] << input_var.getSep() << "sebeta_" << phd.model_terms[i];
- }
- if(interactio != 0) outfile << input_var.getSep() << "beta_SNP_" << phd.model_terms[interaction];
- }
- if (input_var.getNgpreds()==2)
- {
- outfile << input_var.getSep() << "beta_SNP_A1A2"
- << input_var.getSep() << "beta_SNP_A1A1"
- << input_var.getSep() << "sebeta_SNP_A1A2"
- << input_var.getSep() << "sebeta_SNP_a1A1"
- << input_var.getSep() << "chi2_SNP_2df"
- << input_var.getSep() << "beta_SNP_addA1"
- << input_var.getSep() << "sebeta_SNP_addA1"
- << input_var.getSep() << "chi2_SNP_addA1"
- << input_var.getSep() << "beta_SNP_domA1"
- << input_var.getSep() << "sebeta_SNP_domA1"
- << input_var.getSep() << "chi2_SNP_domA1"
- << input_var.getSep() << "beta_SNP_recA1"
- << input_var.getSep() << "sebeta_SNP_recA1"
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/genabel -r 1230
More information about the Genabel-commits
mailing list