aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--Makefile2
-rw-r--r--Makefile.osx13
-rw-r--r--src/bslmmdap.cpp3
-rw-r--r--src/io.cpp31
-rw-r--r--src/io.h2
-rw-r--r--src/lmm.cpp24
-rw-r--r--src/param.h1
-rw-r--r--src/vc.cpp27
8 files changed, 20 insertions, 83 deletions
diff --git a/Makefile b/Makefile
index daa7b21..1b5596a 100644
--- a/Makefile
+++ b/Makefile
@@ -62,7 +62,7 @@ ifdef FORCE_FLOAT
SOURCES += $(SRC_DIR)/param_float.cpp $(SRC_DIR)/gemma_float.cpp $(SRC_DIR)/io_float.cpp $(SRC_DIR)/lm_float.cpp $(SRC_DIR)/vc_float.cpp $(SRC_DIR)/lmm_float.cpp $(SRC_DIR)/mvlmm_float.cpp $(SRC_DIR)/bslmm_float.cpp $(SRC_DIR)/prdt_float.cpp $(SRC_DIR)/mathfunc_float.cpp $(SRC_DIR)/gzstream.cpp $(SRC_DIR)/eigenlib.cpp
HDR += $(SRC_DIR)/param_float.h $(SRC_DIR)/gemma_float.h $(SRC_DIR)/io_float.h $(SRC_DIR)/lm_float.h $(SRC_DIR)/lmm_float.h $(SRC_DIR)/vc_float.h $(SRC_DIR)/mvlmm_float.h $(SRC_DIR)/bslmm_float.h $(SRC_DIR)/prdt_float.h $(SRC_DIR)/mathfunc_float.h $(SRC_DIR)/gzstream.h $(SRC_DIR)/eigenlib.h
else
- SOURCES += $(SRC_DIR)/param.cpp $(SRC_DIR)/gemma.cpp $(SRC_DIR)/io.cpp $(SRC_DIR)/lm.cpp $(SRC_DIR)/lmm.cpp $(SRC_DIR)/vc.cpp $(SRC_DIR)/mvlmm.cpp $(SRC_DIR)/bslmm.cpp $(SRC_DIR)/prdt.cpp $(SRC_DIR)/mathfunc.cpp $(SRC_DIR)/gzstream.cpp $(SRC_DIR)/eigenlib.cpp
+ SOURCES += $(SRC_DIR)/param.cpp $(SRC_DIR)/gemma.cpp $(SRC_DIR)/io.cpp $(SRC_DIR)/lm.cpp $(SRC_DIR)/lmm.cpp $(SRC_DIR)/vc.cpp $(SRC_DIR)/mvlmm.cpp $(SRC_DIR)/bslmm.cpp $(SRC_DIR)/prdt.cpp $(SRC_DIR)/mathfunc.cpp $(SRC_DIR)/gzstream.cpp $(SRC_DIR)/eigenlib.cpp $(SRC_DIR)/ldr.cpp $(SRC_DIR)/bslmmdap.cpp $(SRC_DIR)/logistic.cpp $(SRC_DIR)/varcov.cpp
HDR += $(SRC_DIR)/param.h $(SRC_DIR)/gemma.h $(SRC_DIR)/io.h $(SRC_DIR)/lm.h $(SRC_DIR)/lmm.h $(SRC_DIR)/vc.h $(SRC_DIR)/mvlmm.h $(SRC_DIR)/bslmm.h $(SRC_DIR)/prdt.h $(SRC_DIR)/mathfunc.h $(SRC_DIR)/gzstream.h $(SRC_DIR)/eigenlib.h
endif
diff --git a/Makefile.osx b/Makefile.osx
index e5d5a26..efa2e4c 100644
--- a/Makefile.osx
+++ b/Makefile.osx
@@ -32,7 +32,8 @@ CPP = g++
CPPFLAGS = -F /System/Library/Frameworks -Wall -O3 \
-I/usr/local/opt/gsl@1/include
-LIBS = -lgsl -lgslcblas -lz -L/usr/local/opt/gsl@1/lib
+LIBS = /usr/local/opt/gsl@1/lib/libgsl.a \
+ /usr/local/opt/gsl@1/lib/libgslcblas.a -lz
OUTPUT = $(BIN_DIR)/gemma
@@ -57,8 +58,8 @@ ifdef FORCE_FLOAT
SOURCES += $(SRC_DIR)/param_float.cpp $(SRC_DIR)/gemma_float.cpp $(SRC_DIR)/io_float.cpp $(SRC_DIR)/lm_float.cpp $(SRC_DIR)/vc_float.cpp $(SRC_DIR)/lmm_float.cpp $(SRC_DIR)/mvlmm_float.cpp $(SRC_DIR)/bslmm_float.cpp $(SRC_DIR)/prdt_float.cpp $(SRC_DIR)/mathfunc_float.cpp $(SRC_DIR)/gzstream.cpp $(SRC_DIR)/eigenlib.cpp
HDR += $(SRC_DIR)/param_float.h $(SRC_DIR)/gemma_float.h $(SRC_DIR)/io_float.h $(SRC_DIR)/lm_float.h $(SRC_DIR)/lmm_float.h $(SRC_DIR)/vc_float.h $(SRC_DIR)/mvlmm_float.h $(SRC_DIR)/bslmm_float.h $(SRC_DIR)/prdt_float.h $(SRC_DIR)/mathfunc_float.h $(SRC_DIR)/gzstream.h $(SRC_DIR)/eigenlib.h
else
- SOURCES += $(SRC_DIR)/param.cpp $(SRC_DIR)/gemma.cpp $(SRC_DIR)/io.cpp $(SRC_DIR)/lm.cpp $(SRC_DIR)/lmm.cpp $(SRC_DIR)/vc.cpp $(SRC_DIR)/mvlmm.cpp $(SRC_DIR)/bslmm.cpp $(SRC_DIR)/prdt.cpp $(SRC_DIR)/mathfunc.cpp $(SRC_DIR)/gzstream.cpp $(SRC_DIR)/eigenlib.cpp
- # HDR += $(SRC_DIR)/param.h $(SRC_DIR)/gemma.h $(SRC_DIR)/io.h $(SRC_DIR)/lm.h $(SRC_DIR)/lmm.h $(SRC_DIR)/vc.h $(SRC_DIR)/mvlmm.h $(SRC_DIR)/bslmm.h $(SRC_DIR)/prdt.h $(SRC_DIR)/mathfunc.h $(SRC_DIR)/gzstream.h $(SRC_DIR)/eigenlib.h
+ SOURCES += $(SRC_DIR)/param.cpp $(SRC_DIR)/gemma.cpp $(SRC_DIR)/io.cpp $(SRC_DIR)/lm.cpp $(SRC_DIR)/lmm.cpp $(SRC_DIR)/vc.cpp $(SRC_DIR)/mvlmm.cpp $(SRC_DIR)/bslmm.cpp $(SRC_DIR)/prdt.cpp $(SRC_DIR)/mathfunc.cpp $(SRC_DIR)/gzstream.cpp $(SRC_DIR)/eigenlib.cpp $(SRC_DIR)/ldr.cpp $(SRC_DIR)/bslmmdap.cpp $(SRC_DIR)/logistic.cpp $(SRC_DIR)/varcov.cpp
+ HDR += $(SRC_DIR)/param.h $(SRC_DIR)/gemma.h $(SRC_DIR)/io.h $(SRC_DIR)/lm.h $(SRC_DIR)/lmm.h $(SRC_DIR)/vc.h $(SRC_DIR)/mvlmm.h $(SRC_DIR)/bslmm.h $(SRC_DIR)/prdt.h $(SRC_DIR)/mathfunc.h $(SRC_DIR)/gzstream.h $(SRC_DIR)/eigenlib.h
endif
ifdef WITH_LAPACK
@@ -83,12 +84,6 @@ else
CPPFLAGS += -m64
endif
-ifdef FORCE_DYNAMIC
-else
- CPPFLAGS += -static
-endif
-
-
# all
OBJS = $(SOURCES:.cpp=.o)
diff --git a/src/bslmmdap.cpp b/src/bslmmdap.cpp
index 0bf0e7b..f6e9e9c 100644
--- a/src/bslmmdap.cpp
+++ b/src/bslmmdap.cpp
@@ -42,6 +42,7 @@
#include "logistic.h"
#include "lapack.h"
+#include "io.h"
#ifdef FORCE_FLOAT
#include "param_float.h"
@@ -221,7 +222,7 @@ void ReadFile_cat (const string &file_cat, const vector<string> &vec_rs, gsl_mat
//read header
HEADER header;
!safeGetline(infile, line).eof();
- ReadHeader (line, header);
+ ReadHeader_io (line, header);
//use the header to determine the number of categories
kc=header.catc_col.size(); kd=header.catd_col.size();
diff --git a/src/io.cpp b/src/io.cpp
index 97b29b0..4da1590 100644
--- a/src/io.cpp
+++ b/src/io.cpp
@@ -50,8 +50,6 @@
using namespace std;
-
-
//Print process bar
void ProgressBar (string str, double p, double total)
{
@@ -181,7 +179,7 @@ bool ReadFile_snps_header (const string &file_snps, set<string> &setSnps)
//read header
HEADER header;
!safeGetline(infile, line).eof();
- ReadHeader (line, header);
+ ReadHeader_io (line, header);
if (header.rs_col==0 && (header.chr_col==0 || header.pos_col==0) ) {
cout<<"missing rs id in the hearder"<<endl;
@@ -2501,7 +2499,7 @@ bool bgenKin (const string &file_oxford, vector<int> &indicator_snp, const int k
//read header to determine which column contains which item
-bool ReadHeader (const string &line, HEADER &header)
+bool ReadHeader_io (const string &line, HEADER &header)
{
string rs_ptr[]={"rs","RS","snp","SNP","snps","SNPS","snpid","SNPID","rsid","RSID","MarkerName"};
set<string> rs_set(rs_ptr, rs_ptr+11);
@@ -2645,7 +2643,7 @@ bool ReadFile_cat (const string &file_cat, map<string, size_t> &mapRS2cat, size_
//read header
HEADER header;
!safeGetline(infile, line).eof();
- ReadHeader (line, header);
+ ReadHeader_io (line, header);
//use the header to count the number of categories
n_vc=header.coln;
@@ -2746,7 +2744,7 @@ bool ReadFile_catc (const string &file_cat, map<string, vector<double> > &mapRS2
//read header
HEADER header;
!safeGetline(infile, line).eof();
- ReadHeader (line, header);
+ ReadHeader_io (line, header);
//use the header to count the number of categories
n_cat=header.coln;
@@ -3330,7 +3328,7 @@ bool ReadFile_wsnp (const string &file_wcat, const size_t n_vc, map<string, vect
//read header
HEADER header;
!safeGetline(infile, line).eof();
- ReadHeader (line, header);
+ ReadHeader_io (line, header);
while (!safeGetline(infile, line).eof()) {
if (isBlankLine(line)) {continue;}
@@ -3403,7 +3401,7 @@ void ReadFile_beta (const string &file_beta, const map<string, size_t> &mapRS2ca
//read header
HEADER header;
!safeGetline(infile, line).eof();
- ReadHeader (line, header);
+ ReadHeader_io (line, header);
if (header.n_col==0 ) {
if ( (header.nobs_col==0 && header.nmis_col==0) && (header.ncase_col==0 && header.ncontrol_col==0) ) {
@@ -3534,7 +3532,7 @@ void ReadFile_beta (const string &file_beta, const map<string, double> &mapRS2wA
//read header
HEADER header;
!safeGetline(infile, line).eof();
- ReadHeader (line, header);
+ ReadHeader_io (line, header);
if (header.n_col==0 ) {
if ( (header.nobs_col==0 && header.nmis_col==0) && (header.ncase_col==0 && header.ncontrol_col==0) ) {
@@ -3966,21 +3964,6 @@ void ReadFile_mstudy (const string &file_mstudy, gsl_matrix *Vq_mat, gsl_vector
return;
}
-
-//copied from lmm.cpp; is used in the following function compKtoV
-//map a number 1-(n_cvt+2) to an index between 0 and [(n_c+2)^2+(n_c+2)]/2-1
-size_t GetabIndex (const size_t a, const size_t b, const size_t n_cvt) {
- if (a>n_cvt+2 || b>n_cvt+2 || a<=0 || b<=0) {cout<<"error in GetabIndex."<<endl; return 0;}
- size_t index;
- size_t l, h;
- if (b>a) {l=a; h=b;} else {l=b; h=a;}
-
- size_t n=n_cvt+2;
- index=(2*n-l+2)*(l-1)/2+h-l;
-
- return index;
-}
-
//read reference file
void ReadFile_mref (const string &file_mref, gsl_matrix *S_mat, gsl_matrix *Svar_mat, gsl_vector *s_vec, size_t &ni)
{
diff --git a/src/io.h b/src/io.h
index c1b762d..073812c 100644
--- a/src/io.h
+++ b/src/io.h
@@ -79,7 +79,7 @@ bool CountFileLines (const string &file_input, size_t &n_lines);
bool ReadFile_gene (const string &file_gene, vector<double> &vec_read, vector<SNPINFO> &snpInfo, size_t &ng_total);
-bool ReadHeader (const string &line, HEADER &header);
+bool ReadHeader_io (const string &line, HEADER &header);
bool ReadFile_cat (const string &file_cat, map<string, size_t> &mapRS2cat, size_t &n_vc);
bool ReadFile_mcat (const string &file_mcat, map<string, size_t> &mapRS2cat, size_t &n_vc);
diff --git a/src/lmm.cpp b/src/lmm.cpp
index 044f33c..a707534 100644
--- a/src/lmm.cpp
+++ b/src/lmm.cpp
@@ -183,30 +183,6 @@ void LMM::WriteFiles ()
return;
}
-
-
-
-
-
-
-
-
-
-
-//map a number 1-(n_cvt+2) to an index between 0 and [(n_c+2)^2+(n_c+2)]/2-1
-size_t GetabIndex (const size_t a, const size_t b, const size_t n_cvt) {
- if (a>n_cvt+2 || b>n_cvt+2 || a<=0 || b<=0) {cout<<"error in GetabIndex."<<endl; return 0;}
- size_t index;
- size_t l, h;
- if (b>a) {l=a; h=b;} else {l=b; h=a;}
-
- size_t n=n_cvt+2;
- index=(2*n-l+2)*(l-1)/2+h-l;
-
- return index;
-}
-
-
void CalcPab (const size_t n_cvt, const size_t e_mode, const gsl_vector *Hi_eval, const gsl_matrix *Uab, const gsl_vector *ab, gsl_matrix *Pab)
{
size_t index_ab, index_aw, index_bw, index_ww;
diff --git a/src/param.h b/src/param.h
index 6cb0d97..72b7d85 100644
--- a/src/param.h
+++ b/src/param.h
@@ -307,6 +307,7 @@ public:
void UpdateSNP (const map<string, double> &mapRS2wA);
};
+size_t GetabIndex (const size_t a, const size_t b, const size_t n_cvt);
#endif
diff --git a/src/vc.cpp b/src/vc.cpp
index c0aa40d..1e0c562 100644
--- a/src/vc.cpp
+++ b/src/vc.cpp
@@ -453,7 +453,7 @@ int LogRL_dev12 (const gsl_vector *log_sigma2, void *params, gsl_vector *dev1, g
//read header to determine which column contains which item
-bool ReadHeader (const string &line, HEADER &header)
+bool ReadHeader_vc (const string &line, HEADER &header)
{
string rs_ptr[]={"rs","RS","snp","SNP","snps","SNPS","snpid","SNPID","rsid","RSID"};
set<string> rs_set(rs_ptr, rs_ptr+10);
@@ -586,7 +586,7 @@ void ReadFile_cor (const string &file_cor, const set<string> &setSnps, vector<st
//header
!safeGetline(infile, line).eof();
- ReadHeader (line, header);
+ ReadHeader_vc (line, header);
if (header.n_col==0 ) {
if (header.nobs_col==0 && header.nmis_col==0) {
@@ -700,7 +700,7 @@ void ReadFile_beta (const bool flag_priorscale, const string &file_beta, const m
//read header
HEADER header;
!safeGetline(infile, line).eof();
- ReadHeader (line, header);
+ ReadHeader_vc (line, header);
if (header.n_col==0 ) {
if (header.nobs_col==0 && header.nmis_col==0) {
@@ -869,7 +869,7 @@ void ReadFile_cor (const string &file_cor, const vector<string> &vec_rs, const v
HEADER header;
!safeGetline(infile, line).eof();
- ReadHeader (line, header);
+ ReadHeader_vc (line, header);
while (!safeGetline(infile, line).eof()) {
//do not read cor values this time; upto col_n-1
@@ -1059,25 +1059,6 @@ void ReadFile_cor (const string &file_cor, const vector<string> &vec_rs, const v
return;
}
-
-
-
-
-//copied from lmm.cpp; is used in the following function VCss
-//map a number 1-(n_cvt+2) to an index between 0 and [(n_c+2)^2+(n_c+2)]/2-1
-size_t GetabIndex (const size_t a, const size_t b, const size_t n_cvt) {
- if (a>n_cvt+2 || b>n_cvt+2 || a<=0 || b<=0) {cout<<"error in GetabIndex."<<endl; return 0;}
- size_t index;
- size_t l, h;
- if (b>a) {l=a; h=b;} else {l=b; h=a;}
-
- size_t n=n_cvt+2;
- index=(2*n-l+2)*(l-1)/2+h-l;
-
- return index;
-}
-
-
//use the new method to calculate variance components with summary statistics
//first, use a function CalcS to compute S matrix (where the diagonal elements are part of V(q) ), and then use bootstrap to compute the variance for S, use a set of genotypes, phenotypes, and individual ids, and snp category label
void CalcVCss(const gsl_matrix *Vq, const gsl_matrix *S_mat, const gsl_matrix *Svar_mat, const gsl_vector *q_vec, const gsl_vector *s_vec, const double df, vector<double> &v_pve, vector<double> &v_se_pve, double &pve_total, double &se_pve_total, vector<double> &v_sigma2, vector<double> &v_se_sigma2, vector<double> &v_enrich, vector<double> &v_se_enrich) {