about summary refs log tree commit diff
diff options
context:
space:
mode:
-rw-r--r--src/bslmm.cpp2
-rw-r--r--src/bslmm.h5
-rw-r--r--src/bslmmdap.cpp4
-rw-r--r--src/bslmmdap.h4
-rw-r--r--src/eigenlib.h2
-rw-r--r--src/gzstream.h3
-rw-r--r--src/lapack.h2
-rw-r--r--src/ldr.cpp207
-rw-r--r--src/ldr.h50
-rw-r--r--src/main.cpp24
-rw-r--r--src/varcov.cpp6
-rw-r--r--src/varcov.h6
12 files changed, 59 insertions, 256 deletions
diff --git a/src/bslmm.cpp b/src/bslmm.cpp
index 0bd4234..d295fd8 100644
--- a/src/bslmm.cpp
+++ b/src/bslmm.cpp
@@ -1,6 +1,6 @@
 /*
  Genome-wide Efficient Mixed Model Association (GEMMA)
- Copyright (C) 2011-2017 Xiang Zhou
+ Copyright (C) 2011-2017, Xiang Zhou
  
  This program is free software: you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
diff --git a/src/bslmm.h b/src/bslmm.h
index 863a22d..07aac67 100644
--- a/src/bslmm.h
+++ b/src/bslmm.h
@@ -1,6 +1,6 @@
 /*
  Genome-wide Efficient Mixed Model Association (GEMMA)
- Copyright (C) 2011-2017 Xiang Zhou
+ Copyright (C) 2011-2017, Xiang Zhou
  
  This program is free software: you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
@@ -14,8 +14,7 @@
  
  You should have received a copy of the GNU General Public License
  along with this program. If not, see <http://www.gnu.org/licenses/>.
- */
-
+*/
 
 #ifndef __BSLMM_H__                
 #define __BSLMM_H__
diff --git a/src/bslmmdap.cpp b/src/bslmmdap.cpp
index c66189b..e1b20d8 100644
--- a/src/bslmmdap.cpp
+++ b/src/bslmmdap.cpp
@@ -1,6 +1,6 @@
 /*
  Genome-wide Efficient Mixed Model Association (GEMMA)
- Copyright (C) 2011-2017 Xiang Zhou
+ Copyright (C) 2011-2017, Xiang Zhou
 
  This program is free software: you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
@@ -14,7 +14,7 @@
 
  You should have received a copy of the GNU General Public License
  along with this program. If not, see <http://www.gnu.org/licenses/>.
- */
+*/
 
 #include <iostream>
 #include <fstream>
diff --git a/src/bslmmdap.h b/src/bslmmdap.h
index a9e2d51..7d95db7 100644
--- a/src/bslmmdap.h
+++ b/src/bslmmdap.h
@@ -1,6 +1,6 @@
 /*
  Genome-wide Efficient Mixed Model Association (GEMMA)
- Copyright (C) 2011-2017 Xiang Zhou
+ Copyright (C) 2011-2017, Xiang Zhou
 
  This program is free software: you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
@@ -14,7 +14,7 @@
 
  You should have received a copy of the GNU General Public License
  along with this program. If not, see <http://www.gnu.org/licenses/>.
- */
+*/
 
 #ifndef __BSLMMDAP_H__
 #define __BSLMMDAP_H__
diff --git a/src/eigenlib.h b/src/eigenlib.h
index 493907c..f869786 100644
--- a/src/eigenlib.h
+++ b/src/eigenlib.h
@@ -1,6 +1,6 @@
 /*
     Genome-wide Efficient Mixed Model Association (GEMMA)
-    Copyright (C) 2011-2017 Xiang Zhou
+    Copyright (C) 2011-2017, Xiang Zhou
 
     This program is free software: you can redistribute it and/or modify
     it under the terms of the GNU General Public License as published by
diff --git a/src/gzstream.h b/src/gzstream.h
index 861653f..4e86ad9 100644
--- a/src/gzstream.h
+++ b/src/gzstream.h
@@ -29,7 +29,7 @@
 #ifndef GZSTREAM_H
 #define GZSTREAM_H 1
 
-// standard C++ with new header file names and std:: namespace
+// Standard C++ with new header file names and std::namespace.
 #include <iostream>
 #include <fstream>
 #include <zlib.h>
@@ -118,4 +118,3 @@ public:
 #endif // GZSTREAM_H
 // ============================================================================
 // EOF //
-
diff --git a/src/lapack.h b/src/lapack.h
index 5277b2f..88fc509 100644
--- a/src/lapack.h
+++ b/src/lapack.h
@@ -1,6 +1,6 @@
 /*
     Genome-wide Efficient Mixed Model Association (GEMMA)
-    Copyright (C) 2011-2017 Xiang Zhou
+    Copyright (C) 2011-2017, Xiang Zhou
 
     This program is free software: you can redistribute it and/or modify
     it under the terms of the GNU General Public License as published by
diff --git a/src/ldr.cpp b/src/ldr.cpp
index e28490a..a1e5791 100644
--- a/src/ldr.cpp
+++ b/src/ldr.cpp
@@ -1,6 +1,6 @@
 /*
  Genome-wide Efficient Mixed Model Association (GEMMA)
- Copyright (C) 2011  Xiang Zhou
+ Copyright (C) 2011-2017, Xiang Zhou
  
  This program is free software: you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
@@ -13,8 +13,8 @@
  GNU General Public License for more details.
  
  You should have received a copy of the GNU General Public License
- along with this program.  If not, see <http://www.gnu.org/licenses/>.
- */
+ along with this program. If not, see <http://www.gnu.org/licenses/>.
+*/
 
 #include <iostream>
 #include <fstream>
@@ -39,29 +39,16 @@
 #include "gsl/gsl_roots.h"
 #include "Eigen/Dense"
 
-
-
 #include "lapack.h"
-
-#ifdef FORCE_FLOAT
-#include "param_float.h"
-#include "ldr_float.h"
-#include "lm_float.h"
-#include "mathfunc_float.h"  //for function CenterVector
-#else
 #include "param.h"
 #include "ldr.h"
 #include "lm.h"
 #include "mathfunc.h"
-#endif
 
 using namespace std;
 using namespace Eigen;
 
-
-
-void LDR::CopyFromParam (PARAM &cPar) 
-{
+void LDR::CopyFromParam (PARAM &cPar) {
 	a_mode=cPar.a_mode;
 	d_pace=cPar.d_pace;
 	
@@ -83,176 +70,15 @@ void LDR::CopyFromParam (PARAM &cPar)
 	return;
 }
 
-
-void LDR::CopyToParam (PARAM &cPar) 
-{
-	//cPar.pheno_mean=pheno_mean;
-	//cPar.randseed=randseed;
-	
-	return;
-}
-
-
-/*
-void BSLMM::WriteBV (const gsl_vector *bv) 
-{
-	string file_str;
-	file_str=path_out+"/"+file_out;
-	file_str+=".bv.txt";
-
-	ofstream outfile (file_str.c_str(), ofstream::out);
-	if (!outfile) {cout<<"error writing file: "<<file_str.c_str()<<endl; return;}
-	
-	size_t t=0;
-	for (size_t i=0; i<ni_total; ++i) {
-		if (indicator_idv[i]==0) {
-			outfile<<"NA"<<endl;
-		}		
-		else {
-			outfile<<scientific<<setprecision(6)<<gsl_vector_get(bv, t)<<endl;
-			t++;
-		}
-	}		
-	
-	outfile.clear();	
-	outfile.close();	
+void LDR::CopyToParam (PARAM &cPar) {
 	return;
 }
 
-
-
-
-void BSLMM::WriteParam (vector<pair<double, double> > &beta_g, const gsl_vector *alpha, const size_t w) 
-{
-	string file_str;
-	file_str=path_out+"/"+file_out;
-	file_str+=".param.txt";
-
-	ofstream outfile (file_str.c_str(), ofstream::out);
-	if (!outfile) {cout<<"error writing file: "<<file_str.c_str()<<endl; return;}
-	
-	outfile<<"chr"<<"\t"<<"rs"<<"\t"
-			<<"ps"<<"\t"<<"n_miss"<<"\t"<<"alpha"<<"\t"
-			<<"beta"<<"\t"<<"gamma"<<endl;
-	
-	size_t t=0;
-	for (size_t i=0; i<ns_total; ++i) {
-		if (indicator_snp[i]==0) {continue;}		
-		
-		outfile<<snpInfo[i].chr<<"\t"<<snpInfo[i].rs_number<<"\t"
-		<<snpInfo[i].base_position<<"\t"<<snpInfo[i].n_miss<<"\t";	
-				
-		outfile<<scientific<<setprecision(6)<<gsl_vector_get(alpha, t)<<"\t";
-		if (beta_g[t].second!=0) {
-			outfile<<beta_g[t].first/beta_g[t].second<<"\t"<<beta_g[t].second/(double)w<<endl;
-		}
-		else {
-			outfile<<0.0<<"\t"<<0.0<<endl;
-		}
-		t++;
-	}		
-	
-	outfile.clear();	
-	outfile.close();	
-	return;
-}
-
-
-void BSLMM::WriteParam (const gsl_vector *alpha) 
-{
-	string file_str;
-	file_str=path_out+"/"+file_out;
-	file_str+=".param.txt";
-
-	ofstream outfile (file_str.c_str(), ofstream::out);
-	if (!outfile) {cout<<"error writing file: "<<file_str.c_str()<<endl; return;}
-	
-	outfile<<"chr"<<"\t"<<"rs"<<"\t"
-			<<"ps"<<"\t"<<"n_miss"<<"\t"<<"alpha"<<"\t"
-			<<"beta"<<"\t"<<"gamma"<<endl;
-	
-	size_t t=0;
-	for (size_t i=0; i<ns_total; ++i) {
-		if (indicator_snp[i]==0) {continue;}		
-
-		outfile<<snpInfo[i].chr<<"\t"<<snpInfo[i].rs_number<<"\t"
-				<<snpInfo[i].base_position<<"\t"<<snpInfo[i].n_miss<<"\t";				
-		outfile<<scientific<<setprecision(6)<<gsl_vector_get(alpha, t)<<"\t";
-		outfile<<0.0<<"\t"<<0.0<<endl;
-		t++;
-	}		
-	
-	outfile.clear();	
-	outfile.close();	
-	return;
-}
-
-
-void BSLMM::WriteResult (const int flag, const gsl_matrix *Result_hyp, const gsl_matrix *Result_gamma, const size_t w_col) 
-{
-	string file_gamma, file_hyp;
-	file_gamma=path_out+"/"+file_out;
-	file_gamma+=".gamma.txt";
-	file_hyp=path_out+"/"+file_out;
-	file_hyp+=".hyp.txt";
-
-	ofstream outfile_gamma, outfile_hyp;
-		
-	if (flag==0) {
-		outfile_gamma.open (file_gamma.c_str(), ofstream::out);
-		outfile_hyp.open (file_hyp.c_str(), ofstream::out);
-		if (!outfile_gamma) {cout<<"error writing file: "<<file_gamma<<endl; return;}
-		if (!outfile_hyp) {cout<<"error writing file: "<<file_hyp<<endl; return;}
-		
-		outfile_hyp<<"h \t pve \t rho \t pge \t pi \t n_gamma"<<endl;
-		
-		for (size_t i=0; i<s_max; ++i) {
-			outfile_gamma<<"s"<<i<<"\t";
-		}
-		outfile_gamma<<endl;
-	}
-	else {
-		outfile_gamma.open (file_gamma.c_str(), ofstream::app);
-		outfile_hyp.open (file_hyp.c_str(), ofstream::app);
-		if (!outfile_gamma) {cout<<"error writing file: "<<file_gamma<<endl; return;}
-		if (!outfile_hyp) {cout<<"error writing file: "<<file_hyp<<endl; return;}
-		
-		size_t w;
-		if (w_col==0) {w=w_pace;}
-		else {w=w_col;}
-		
-		for (size_t i=0; i<w; ++i) {
-			outfile_hyp<<scientific;
-			for (size_t j=0; j<4; ++j) {
-				outfile_hyp<<setprecision(6)<<gsl_matrix_get (Result_hyp, i, j)<<"\t";
-			}
-			outfile_hyp<<setprecision(6)<<exp(gsl_matrix_get (Result_hyp, i, 4))<<"\t";
-			outfile_hyp<<(int)gsl_matrix_get (Result_hyp, i, 5)<<"\t";
-			outfile_hyp<<endl;
-		}
-		
-		for (size_t i=0; i<w; ++i) {
-			for (size_t j=0; j<s_max; ++j) {
-				outfile_gamma<<(int)gsl_matrix_get (Result_gamma, i, j)<<"\t";
-			}
-			outfile_gamma<<endl;
-		}
-		
-	}
-	
-	outfile_hyp.close();
-	outfile_hyp.clear();
-	outfile_gamma.close();
-	outfile_gamma.clear();	
-	return;
-}
-
-*/
-
-//X is a p by n matrix
-void LDR::VB (const vector<vector<unsigned char> > &Xt, const gsl_matrix *W_gsl, const gsl_vector *y_gsl)
-{ 
-  //save gsl_vector and gsl_matrix into eigen library formats
+//X is a p by n matrix.
+void LDR::VB (const vector<vector<unsigned char> > &Xt,
+	      const gsl_matrix *W_gsl, const gsl_vector *y_gsl) {
+  
+  // Save gsl_vector and gsl_matrix into Eigen library formats.
   MatrixXd W(W_gsl->size1, W_gsl->size2);
   VectorXd y(y_gsl->size);
   VectorXd x_col(y_gsl->size);
@@ -266,7 +92,7 @@ void LDR::VB (const vector<vector<unsigned char> > &Xt, const gsl_matrix *W_gsl,
     }
   }
 
-  //initial VB values by lm
+  // Initial VB values by lm.
   cout<<indicator_snp[0]<<" "<<indicator_snp[1]<<" "<<indicator_snp[2]<<endl;
   uchar_matrix_get_row (Xt, 0, x_col);
 
@@ -274,12 +100,11 @@ void LDR::VB (const vector<vector<unsigned char> > &Xt, const gsl_matrix *W_gsl,
     cout<<x_col(j)<<endl;
   }
 
+  // Run VB iterations.
+  // TO DO.
 
-  //run VB iterations
-
-
-
-  //save results
-
+  // Save results.
+  // TO DO.
+  
   return;
 }
diff --git a/src/ldr.h b/src/ldr.h
index d81048c..ceb08cf 100644
--- a/src/ldr.h
+++ b/src/ldr.h
@@ -1,6 +1,6 @@
 /*
  Genome-wide Efficient Mixed Model Association (GEMMA)
- Copyright (C) 2011  Xiang Zhou
+ Copyright (C) 2011-2017, Xiang Zhou
  
  This program is free software: you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
@@ -13,9 +13,8 @@
  GNU General Public License for more details.
  
  You should have received a copy of the GNU General Public License
- along with this program.  If not, see <http://www.gnu.org/licenses/>.
- */
-
+ along with this program. If not, see <http://www.gnu.org/licenses/>.
+*/
 
 #ifndef __LDR_H__                
 #define __LDR_H__
@@ -24,25 +23,14 @@
 #include <map>
 #include <gsl/gsl_rng.h>
 #include <gsl/gsl_randist.h>
-
-#ifdef FORCE_FLOAT
-#include "param_float.h"
-#else
 #include "param.h"
-#endif
-
 
 using namespace std;
 
-
-
-
-
-
 class LDR {
 
 public:	
-	// IO related parameters
+	// IO-related parameters.
 	int a_mode;	
 	size_t d_pace;
 	
@@ -51,27 +39,33 @@ public:
 	string file_out;
 	string path_out;
 		
-	// Summary statistics
-	size_t ni_total, ns_total;	//number of total individuals and snps
-	size_t ni_test, ns_test;	//number of individuals and snps used for analysis
-	size_t n_cvt;				//number of covariates
-	vector<int> indicator_idv;				//indicator for individuals (phenotypes), 0 missing, 1 available for analysis
-	vector<int> indicator_snp;				//sequence indicator for SNPs: 0 ignored because of (a) maf, (b) miss, (c) non-poly; 1 available for analysis
+	// Summary statistics.
+	size_t ni_total, ns_total; // Total number of individuals & SNPs.
+	size_t ni_test, ns_test;   // Number of individuals & SNPs used
+                                   // for analysis
+	size_t n_cvt;		   // Number of covariates.
+
+        // Indicator for individuals (phenotypes): 0 missing, 1
+        // available for analysis.
+	vector<int> indicator_idv;
+
+       	// Sequence indicator for SNPs: 0 ignored because of (a) maf,
+       	// (b) miss, (c) non-poly; 1 available for analysis.
+	vector<int> indicator_snp;
 	
-	vector<SNPINFO> snpInfo;		//record SNP information
+	vector<SNPINFO> snpInfo; // Record SNP information.
 	
-	// Not included in PARAM
+	// Not included in PARAM.
 	gsl_rng *gsl_r; 	
 	
-	// Main Functions
+	// Main functions.
 	void CopyFromParam (PARAM &cPar);
 	void CopyToParam (PARAM &cPar);
 	
-	void VB(const vector<vector<unsigned char> > &Xt, const gsl_matrix *W_gsl, const gsl_vector *y_gsl);
+	void VB(const vector<vector<unsigned char> > &Xt,
+		const gsl_matrix *W_gsl, const gsl_vector *y_gsl);
 };
 
-
-
 #endif
 
 
diff --git a/src/main.cpp b/src/main.cpp
index e1fb336..b7ac884 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -1,6 +1,6 @@
 /*
-	Genome-wide Efficient Mixed Model Association (GEMMA)
-    Copyright (C) 2011  Xiang Zhou
+    Genome-wide Efficient Mixed Model Association (GEMMA)
+    Copyright (C) 2011-2017, Xiang Zhou
 
     This program is free software: you can redistribute it and/or modify
     it under the terms of the GNU General Public License as published by
@@ -13,7 +13,7 @@
     GNU General Public License for more details.
 
     You should have received a copy of the GNU General Public License
-    along with this program.  If not, see <http://www.gnu.org/licenses/>.
+    along with this program. If not, see <http://www.gnu.org/licenses/>.
 */
 
 #include <iostream>
@@ -21,21 +21,11 @@
 #include <sstream>
 #include <sys/stat.h>
 #include <sys/types.h>
-
-#include "param.h"
-
-#ifdef FORCE_FLOAT
-#include "gemma_float.h"
-#else
 #include "gemma.h"
-#endif
 
 using namespace std;
 
-
-
-int main(int argc, char * argv[])
-{ 	
+int main(int argc, char * argv[]) { 	
 	GEMMA cGemma;	
 	PARAM cPar;
 
@@ -79,8 +69,4 @@ int main(int argc, char * argv[])
 	
 	cGemma.WriteLog(argc, argv, cPar);
 	
-    return EXIT_SUCCESS;                                                          
-}
-
-
- 
+    return EXIT_SUCCESS;                                                       }
diff --git a/src/varcov.cpp b/src/varcov.cpp
index 66cf3c4..9d39c6c 100644
--- a/src/varcov.cpp
+++ b/src/varcov.cpp
@@ -1,6 +1,6 @@
 /*
-	Genome-wide Efficient Mixed Model Association (GEMMA)
-    Copyright (C) 2011-2017 Xiang Zhou
+    Genome-wide Efficient Mixed Model Association (GEMMA)
+    Copyright (C) 2011-2017, Xiang Zhou
 
     This program is free software: you can redistribute it and/or modify
     it under the terms of the GNU General Public License as published by
@@ -13,7 +13,7 @@
     GNU General Public License for more details.
 
     You should have received a copy of the GNU General Public License
-    along with this program.  If not, see <http://www.gnu.org/licenses/>.
+    along with this program. If not, see <http://www.gnu.org/licenses/>.
 */
 
 #include <iostream>
diff --git a/src/varcov.h b/src/varcov.h
index a7fcd9c..291f16c 100644
--- a/src/varcov.h
+++ b/src/varcov.h
@@ -1,6 +1,6 @@
 /*
-	Genome-wide Efficient Mixed Model Association (GEMMA)
-    Copyright (C) 2011-2017 Xiang Zhou
+    Genome-wide Efficient Mixed Model Association (GEMMA)
+    Copyright (C) 2011-2017, Xiang Zhou
 
     This program is free software: you can redistribute it and/or modify
     it under the terms of the GNU General Public License as published by
@@ -13,7 +13,7 @@
     GNU General Public License for more details.
 
     You should have received a copy of the GNU General Public License
-    along with this program.  If not, see <http://www.gnu.org/licenses/>.
+    along with this program. If not, see <http://www.gnu.org/licenses/>.
 */
 
 #ifndef __VARCOV_H__