[Genabel-commits] r1433 - pkg/ProbABEL/src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Nov 28 20:48:58 CET 2013
Author: maartenk
Date: 2013-11-28 20:48:57 +0100 (Thu, 28 Nov 2013)
New Revision: 1433
Added:
pkg/ProbABEL/src/main_functions_dump.cpp
pkg/ProbABEL/src/main_functions_dump.h
Modified:
pkg/ProbABEL/src/Makefile.am
pkg/ProbABEL/src/coxph_data.cpp
pkg/ProbABEL/src/coxph_data.h
pkg/ProbABEL/src/main.cpp
pkg/ProbABEL/src/regdata.cpp
pkg/ProbABEL/src/regdata.h
Log:
Modified: pkg/ProbABEL/src/Makefile.am
===================================================================
--- pkg/ProbABEL/src/Makefile.am 2013-11-28 12:59:57 UTC (rev 1432)
+++ pkg/ProbABEL/src/Makefile.am 2013-11-28 19:48:57 UTC (rev 1433)
@@ -7,7 +7,7 @@
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 \
- maskedmatrix.cpp maskedmatrix.h reg1.cpp
+ maskedmatrix.cpp maskedmatrix.h reg1.cpp main_functions_dump.h main_functions_dump.cpp
EIGENFILES = eigen_mematrix.h eigen_mematrix.cpp
Modified: pkg/ProbABEL/src/coxph_data.cpp
===================================================================
--- pkg/ProbABEL/src/coxph_data.cpp 2013-11-28 12:59:57 UTC (rev 1432)
+++ pkg/ProbABEL/src/coxph_data.cpp 2013-11-28 19:48:57 UTC (rev 1433)
@@ -55,6 +55,8 @@
strata = obj.strata;
X = obj.X;
order = obj.order;
+ gcount =0;
+ freq =0;
masked_data = new unsigned short int[nids];
for (int i = 0; i < nids; i++)
@@ -185,43 +187,61 @@
delete[] passed_sorted;
}
-void coxph_data::update_snp(gendata &gend, const int snpnum)
-{
- /**
- * This is the main part of the fix of bug #1846
- * (C) of the fix:
- * UMC St Radboud Nijmegen,
- * Dept of Epidemiology & Biostatistics,
- * led by Prof. B. Kiemeney
- *
- * Note this sorts by "order"!!!
- * Here we deal with transposed X, hence last two arguments are swapped
- * compared to the other 'update_snp'
- * Also, the starting column-1 is not necessary for cox X therefore
- * 'ncov-j' changes to 'ncov-j-1'
- **/
+void coxph_data::update_snp(gendata &gend, const int snpnum) {
+ /**
+ * This is the main part of the fix of bug #1846
+ * (C) of the fix:
+ * UMC St Radboud Nijmegen,
+ * Dept of Epidemiology & Biostatistics,
+ * led by Prof. B. Kiemeney
+ *
+ * Note this sorts by "order"!!!
+ * Here we deal with transposed X, hence last two arguments are swapped
+ * compared to the other 'update_snp'
+ * Also, the starting column-1 is not necessary for cox X therefore
+ * 'ncov-j' changes to 'ncov-j-1'
+ **/
+ //reset counter for frequency since it is a new snp
+ gcount=0;
+ freq=0.0;
+ for (int j = 0; j < ngpreds; j++) {
+ double *snpdata = new double[nids];
+ for (int i = 0; i < nids; i++) {
+ masked_data[i] = 0;
+ }
- for (int j = 0; j < ngpreds; j++)
- {
- double *snpdata = new double[nids];
- for (int i = 0; i < nids; i++)
- {
- masked_data[i] = 0;
- }
+ gend.get_var(snpnum * ngpreds + j, snpdata);
- gend.get_var(snpnum * ngpreds + j, snpdata);
+ for (int i = 0; i < nids; i++) {
+ X.put(snpdata[i], (ncov - j - 1), order[i]);
+ if (std::isnan(snpdata[i])) {
+ masked_data[order[i]] = 1;
+ //snp not masked
+ } else {
+ // checck for first predicor
+ if (j == 0) {
+ gcount++;
+ if (ngpreds == 1) {
+ freq += snpdata[i] * 0.5;
+ } else if (ngpreds == 2) {
+ freq += snpdata[i];
+ }
+ } else if (j == 1) {
+ // add second genotype in two predicor data form
+ freq += snpdata[i] * 0.5;
+ }
+ } //end std::isnan(snpdata[i]) snp
- for (int i = 0; i < nids; i++)
- {
- X.put(snpdata[i], (ncov - j - 1), order[i]);
- if (std::isnan(snpdata[i]))
- masked_data[order[i]] = 1;
- }
- delete[] snpdata;
- }
+ } //end i for loop
+
+ delete[] snpdata;
+ } //end ngpreds loop
+ freq /= static_cast<double>(gcount); // Allele frequency
+
}
+
void coxph_data::remove_snp_from_X()
{
// update_snp() adds SNP information to the design matrix. This
Modified: pkg/ProbABEL/src/coxph_data.h
===================================================================
--- pkg/ProbABEL/src/coxph_data.h 2013-11-28 12:59:57 UTC (rev 1432)
+++ pkg/ProbABEL/src/coxph_data.h 2013-11-28 19:48:57 UTC (rev 1433)
@@ -42,6 +42,8 @@
int nids;
int ncov;
int ngpreds;
+ unsigned int gcount;
+ double freq;
mematrix<double> weights;
mematrix<double> stime;
mematrix<int> sstat;
Modified: pkg/ProbABEL/src/main.cpp
===================================================================
--- pkg/ProbABEL/src/main.cpp 2013-11-28 12:59:57 UTC (rev 1432)
+++ pkg/ProbABEL/src/main.cpp 2013-11-28 19:48:57 UTC (rev 1433)
@@ -54,6 +54,7 @@
#include "reg1.h"
#include "command_line_settings.h"
#include "coxph_data.h"
+#include "main_functions_dump.h"
/// Maximum number of iterations we allow the regression to run
#define MAXITER 20
@@ -64,411 +65,8 @@
-/**
- * Send a progress update (a percentage) to stdout so that the user
- * has a rough indication of the percentage of SNPs that has already
- * been completed.
- *
- * @param csnp Number of the SNP that is currently being analysed.
- * @param nsnps Total number of SNPs
- */
-void update_progress_to_cmd_line(int csnp, int nsnps)
-{
- std::cout << setprecision(2) << fixed;
- if (csnp % 1000 == 0)
- {
- if (csnp == 0)
- {
- std::cout << "Analysis: "
- << 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)
- << "%...";
- }
- std::cout.flush();
- }
- std::cout << setprecision(6);
-}
-
/**
- * Open an output file for each model when using probability data
- * (npgreds == 2). This function creates the _2df.out.txt etc. files.
- *
- * @param outfile Vector of output streams
- * @param outfilename_str Basename of the outputfiles.
- */
-void open_files_for_output(std::vector<std::ofstream*>& outfile,
- std::string& outfilename_str)
-{
- //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" };
-
- for (int i = 0; i < amount_of_files; i++)
- {
- outfile.push_back(new std::ofstream());
- outfile[i]->open((filenames[i]).c_str());
- if (!outfile[i]->is_open())
- {
- std::cerr << "Cannot open file for writing: "
- << filenames[i]
- << "\n";
- exit(1);
- }
- }
-}
-
-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());
-
- 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)
- {
- std::cerr << "error: Interaction parameter is out of range "
- << "(interaction="
- << input_var.getInteraction()
- << ") \n";
- exit(1);
- }
-
- return interaction_cox;
-}
-
-
-/**
- * Load the inverse variance-covariance matrix into an InvSigma object.
- *
- * @param input_var Object containing the values of the various
- * command line options.
- * @param phd Object with phenotype data
- * @param invvarmatrix The object of type masked_matrix in which the
- * inverse variance-covariance matrix is returned.
- */
-void loadInvSigma(cmdvars& input_var, phedata& phd, masked_matrix& invvarmatrix)
-{
- std::cout << "You are running mmscore...\n";
- InvSigma inv(input_var.getInverseFilename(), &phd);
- // invvarmatrix = inv.get_matrix();
- //double par = 1.; //var(phd.Y)*phd.nids/(phd.nids-phd.ncov-1);
- invvarmatrix.set_matrix(inv.get_matrix()); // = invvarmatrix * par;
- std::cout << " loaded InvSigma...\n" << std::flush;
-}
-
-
-/**
- * Create the first part of the output file header.
- *
- * \param outfile Vector of output file streams. Contains the streams
- * of the output file(s). One file when using dosage data (ngpreds==1)
- * and one for each genetic model in case probabilities are used
- * (ngpreds==2).
- * \param input_var Object containing the values of the various
- * command line options.
- * \param phd Object with phenotype data
- */
-void create_start_of_header(std::vector<std::ofstream*>& outfile,
- cmdvars& input_var, phedata& phd)
-{
- 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";
- if (input_var.getChrom() != "-1")
- (*outfile[i]) << input_var.getSep() << "chrom";
- if (input_var.getMapfilename() != NULL)
- (*outfile[i]) << input_var.getSep() << "position";
- }
-
- if (input_var.getAllcov()) //All covariates in output
- {
- 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];
- }
-}
-
-
-/**
- * Create the header of the output file(s).
- *
- * \param outfile vector of output file streams. Contains the streams
- * of the output file(s). One file when using dosage data (ngpreds==1)
- * and one for each genetic model in case probabilities are used
- * (ngpreds==2).
- * \param input_var object containing the values of the various
- * command line options.
- * \param phd object with phenotype data
- * \param interaction_cox are we using the Cox model with interaction?
- */
-void create_header(std::vector<std::ofstream*>& outfile, cmdvars& input_var,
- phedata& phd, int& interaction_cox)
-{
- create_start_of_header(outfile, input_var, phd);
-
- if (input_var.getNgpreds() == 1) // dose data: only additive model
- {
- *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_"
- << phd.model_terms[interaction_cox]
- << input_var.getSep()
- << "sebeta_SNP_"
- << phd.model_terms[interaction_cox];
- }
-
- if (input_var.getInverseFilename() == NULL)
- {
- //Han Chen
-#if !COXPH
- if (input_var.getInteraction() != 0 && !input_var.getAllcov())
- {
- *outfile[0] << input_var.getSep()
- << "cov_SNP_int_SNP_"
- << phd.model_terms[interaction_cox];
- }
-#endif
- }
- *outfile[0] << input_var.getSep() << "chi2_SNP";
- *outfile[0] << "\n";
- } // ngpreds == 1
- else if (input_var.getNgpreds() == 2) // prob data: all models
- {
- *outfile[0] << input_var.getSep()
- << "beta_SNP_A1A2"
- << input_var.getSep()
- << "sebeta_SNP_A1A2"
- << input_var.getSep()
- << "beta_SNP_A1A1"
- << 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()
- << "sebeta_SNP_domA1";
- *outfile[3] << input_var.getSep()
- << "beta_SNP_recA1"
- << input_var.getSep()
- << "sebeta_SNP_recA1";
- *outfile[4] << input_var.getSep()
- << "beta_SNP_odomA1"
- << input_var.getSep()
- << "sebeta_SNP_odomA1";
- 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];
-#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];
- }
-#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];
- //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];
- }
-#endif
- //Oct 26, 2009
- }
- }
- *outfile[0] << input_var.getSep() << "chi2_SNP_2df\n"; // "loglik\n";
- *outfile[1] << input_var.getSep() << "chi2_SNP_A1\n"; // "loglik\n";
- *outfile[2] << input_var.getSep() << "chi2_SNP_domA1\n";// "loglik\n";
- *outfile[3] << input_var.getSep() << "chi2_SNP_recA1\n";// "loglik\n";
- *outfile[4] << input_var.getSep() << "chi2_SNP_odomA1\n"; // "loglik\n";
- } // End: ngpreds == 2
- else
- {
- cerr << "Error: create_header(): ngpreds != 1 or 2.\n";
- }
-}
-
-
-/**
- * Write the information from the mlinfo file to the output file(s).
- *
- * \param outfile Vector of output file(s)
- * \param file index number identifying the file in the vector of files
- * \param mli mlinfo object
- * \param csnp number of the SNP that is currently being analysed
- * \param input_var object containing the information of the options
- * specified on the command line
- * \param gcount The number of non-NaN genotypes
- * \param freq The allele frequency based on the non-NaN genotypes
- */
-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];
- }
-}
-
-
-/**
- * Get the position within a (row or column) vector (the index) where
- * a \f$ \beta \f$ (or \f$ se_{\beta} \f$) starts. This is basically a
- * matter of counting backwards from the end of the vector/list.
- *
- * @param input_var Object containing the values of the various
- * command line options.
- * @param model Number of the genetic model (additive, etc)
- * @param number_of_rows_or_columns Total number of rows or columns in
- * the vector.
- *
- * @return Start position of beta for this model
- */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;
-}
-
-
-/**
* Main routine. The main logic of ProbABEL can be found here
*
* \param argc Number of command line arguments
@@ -632,43 +230,10 @@
for (int csnp = 0; csnp < nsnps; csnp++)
{
rgd.update_snp(gtd, csnp);
- double freq = 0.;
- unsigned int gcount = 0;
- double *snpdata1 = new double[gtd.nids];
- double *snpdata2 = new double[gtd.nids];
- if (input_var.getNgpreds() == 2) // Two predictors (probs)
- {
- //freq = ((gtd.G).column_mean(csnp*2)*2. +
- // (gtd.G).column_mean(csnp*2+1))/2.;
- gtd.get_var(csnp * 2, snpdata1);
- gtd.get_var(csnp * 2 + 1, snpdata2);
- for (unsigned int ii = 0; ii < gtd.nids; ii++)
- {
- if (!std::isnan(snpdata1[ii]) && !std::isnan(snpdata2[ii]))
- {
- gcount++;
- freq += snpdata1[ii] + snpdata2[ii] * 0.5;
- }
- }
- }
- else // Only one predictor (dosage data)
- {
- // freq = (gtd.G).column_mean(csnp)/2.;
- gtd.get_var(csnp, snpdata1);
- for (unsigned int ii = 0; ii < gtd.nids; ii++)
- {
- if (!std::isnan(snpdata1[ii]))
- {
- gcount++;
- freq += snpdata1[ii] * 0.5;
- }
- }
- }
- freq /= static_cast<double>(gcount); // Allele frequency
int poly = 1;
- if (fabs(freq) < 1.e-16 || fabs(1. - freq) < 1.e-16)
+ if (fabs(rgd.freq) < 1.e-16 || fabs(1. - rgd.freq) < 1.e-16)
{
poly = 0;
}
@@ -685,13 +250,13 @@
for (unsigned int file = 0; file < outfile.size(); file++)
{
write_mlinfo(outfile, file, mli, csnp, input_var,
- gcount, freq);
+ rgd.gcount, rgd.freq);
}
} else
{
// Dosage data: only additive model
int file = 0;
- write_mlinfo(outfile, file, mli, csnp, input_var, gcount, freq);
+ write_mlinfo(outfile, file, mli, csnp, input_var, rgd.gcount, rgd.freq);
maxmod = 1; // We can only calculate the additive
// model with dosage data
}
@@ -801,7 +366,7 @@
if (input_var.getScore() == 0)
{
double loglik = rd.loglik;
- if (gcount != gtd.nids)
+ if (rgd.gcount != gtd.nids)
{
// If SNP data is missing we didn't
// correctly compute the null likelihood
@@ -992,8 +557,6 @@
}
update_progress_to_cmd_line(csnp, nsnps);
- delete[] snpdata1;
- delete[] snpdata2;
} // END for loop over all SNPs
Added: pkg/ProbABEL/src/main_functions_dump.cpp
===================================================================
--- pkg/ProbABEL/src/main_functions_dump.cpp (rev 0)
+++ pkg/ProbABEL/src/main_functions_dump.cpp 2013-11-28 19:48:57 UTC (rev 1433)
@@ -0,0 +1,431 @@
+/*
+ * main_functions_dump.cpp
+ *
+ * Created on: Nov 27, 2013
+ * Author: mkooyman
+ */
+
+
+#include <stdio.h>
+#include <iostream>
+#include <cstdlib>
+#include <fstream>
+#include <sstream>
+#include <string>
+#include <iomanip>
+#include <vector>
+
+#include "maskedmatrix.h"
+#include "phedata.h"
+#include "data.h"
+#include "command_line_settings.h"
+
+/**
+ * Send a progress update (a percentage) to stdout so that the user
+ * has a rough indication of the percentage of SNPs that has already
+ * been completed.
+ *
+ * @param csnp Number of the SNP that is currently being analysed.
+ * @param nsnps Total number of SNPs
+ */
+void update_progress_to_cmd_line(int csnp, int nsnps)
+{
+ std::cout << std::setprecision(2) << std::fixed;
+
+ if (csnp % 1000 == 0)
+ {
+ if (csnp == 0)
+ {
+ std::cout << "Analysis: "
+ << std::setw(5)
+ << 100. * static_cast<double>(csnp)
+ / static_cast<double>(nsnps)
+ << "%...";
+ }
+ else
+ {
+ std::cout << "\b\b\b\b\b\b\b\b\b"
+ << std::setw(5)
+ << 100. * static_cast<double>(csnp)
+ / static_cast<double>(nsnps)
+ << "%...";
+ }
+ std::cout.flush();
+ }
+ std::cout << std::setprecision(6);
+}
+
+/**
+ * Open an output file for each model when using probability data
+ * (npgreds == 2). This function creates the _2df.out.txt etc. files.
+ *
+ * @param outfile Vector of output streams
+ * @param outfilename_str Basename of the outputfiles.
+ */
+void open_files_for_output(std::vector<std::ofstream*>& outfile,
+ std::string& outfilename_str)
+{
+ //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" };
+
+ for (int i = 0; i < amount_of_files; i++)
+ {
+ outfile.push_back(new std::ofstream());
+ outfile[i]->open((filenames[i]).c_str());
+ if (!outfile[i]->is_open())
+ {
+ std::cerr << "Cannot open file for writing: "
+ << filenames[i]
+ << "\n";
+ exit(1);
+ }
+ }
+}
+
+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());
+
+ 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)
+ {
+ std::cerr << "error: Interaction parameter is out of range "
+ << "(interaction="
+ << input_var.getInteraction()
+ << ") \n";
+ exit(1);
+ }
+
+ return interaction_cox;
+}
+
+
+/**
+ * Load the inverse variance-covariance matrix into an InvSigma object.
+ *
+ * @param input_var Object containing the values of the various
+ * command line options.
+ * @param phd Object with phenotype data
+ * @param invvarmatrix The object of type masked_matrix in which the
+ * inverse variance-covariance matrix is returned.
+ */
+void loadInvSigma(cmdvars& input_var, phedata& phd, masked_matrix& invvarmatrix)
+{
+ std::cout << "You are running mmscore...\n";
+ InvSigma inv(input_var.getInverseFilename(), &phd);
+ // invvarmatrix = inv.get_matrix();
+ //double par = 1.; //var(phd.Y)*phd.nids/(phd.nids-phd.ncov-1);
+ invvarmatrix.set_matrix(inv.get_matrix()); // = invvarmatrix * par;
+ std::cout << " loaded InvSigma...\n" << std::flush;
+}
+
+
+/**
+ * Create the first part of the output file header.
+ *
+ * \param outfile Vector of output file streams. Contains the streams
+ * of the output file(s). One file when using dosage data (ngpreds==1)
+ * and one for each genetic model in case probabilities are used
+ * (ngpreds==2).
+ * \param input_var Object containing the values of the various
+ * command line options.
+ * \param phd Object with phenotype data
+ */
+void create_start_of_header(std::vector<std::ofstream*>& outfile,
+ cmdvars& input_var, phedata& phd)
+{
+ 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";
+ if (input_var.getChrom() != "-1")
+ (*outfile[i]) << input_var.getSep() << "chrom";
+ if (input_var.getMapfilename() != NULL)
+ (*outfile[i]) << input_var.getSep() << "position";
+ }
+
+ if (input_var.getAllcov()) //All covariates in output
+ {
+ 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];
+ }
+}
+
+
+/**
+ * Create the header of the output file(s).
+ *
+ * \param outfile vector of output file streams. Contains the streams
+ * of the output file(s). One file when using dosage data (ngpreds==1)
+ * and one for each genetic model in case probabilities are used
+ * (ngpreds==2).
+ * \param input_var object containing the values of the various
+ * command line options.
+ * \param phd object with phenotype data
+ * \param interaction_cox are we using the Cox model with interaction?
+ */
+void create_header(std::vector<std::ofstream*>& outfile, cmdvars& input_var,
+ phedata& phd, int& interaction_cox)
+{
+ create_start_of_header(outfile, input_var, phd);
+
+ if (input_var.getNgpreds() == 1) // dose data: only additive model
+ {
+ *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_"
+ << phd.model_terms[interaction_cox]
+ << input_var.getSep()
+ << "sebeta_SNP_"
+ << phd.model_terms[interaction_cox];
+ }
+
+ if (input_var.getInverseFilename() == NULL)
+ {
+ //Han Chen
+#if !COXPH
+ if (input_var.getInteraction() != 0 && !input_var.getAllcov())
+ {
+ *outfile[0] << input_var.getSep()
+ << "cov_SNP_int_SNP_"
+ << phd.model_terms[interaction_cox];
+ }
+#endif
+ }
+ *outfile[0] << input_var.getSep() << "chi2_SNP";
+ *outfile[0] << "\n";
+ } // ngpreds == 1
+ else if (input_var.getNgpreds() == 2) // prob data: all models
+ {
+ *outfile[0] << input_var.getSep()
+ << "beta_SNP_A1A2"
+ << input_var.getSep()
+ << "sebeta_SNP_A1A2"
+ << input_var.getSep()
+ << "beta_SNP_A1A1"
+ << 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()
+ << "sebeta_SNP_domA1";
+ *outfile[3] << input_var.getSep()
+ << "beta_SNP_recA1"
+ << input_var.getSep()
+ << "sebeta_SNP_recA1";
+ *outfile[4] << input_var.getSep()
+ << "beta_SNP_odomA1"
+ << input_var.getSep()
+ << "sebeta_SNP_odomA1";
+ 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];
+#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];
+ }
+#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];
+ //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];
+ }
+#endif
+ //Oct 26, 2009
+ }
+ }
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/genabel -r 1433
More information about the Genabel-commits
mailing list