version 0.95alpha

1 files changed, 146 insertions, 0 deletions

diff --git a/src/bslmm.h b/src/bslmm.h
new file mode 100644
index 0000000..8b5edc7
--- /dev/null
+++ b/src/bslmm.h
@@ -0,0 +1,146 @@
+/*
+ Genome-wide Efficient Mixed Model Association (GEMMA)
+ Copyright (C) 2011  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
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+ 
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ 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/>.
+ */
+
+
+#ifndef __BSLMM_H__                
+#define __BSLMM_H__
+
+#include <vector>
+#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 BSLMM {
+
+public:	
+	// IO related parameters
+	int a_mode;	
+	size_t d_pace;
+	
+	string file_bfile;
+	string file_geno;
+	string file_out;
+	string path_out;
+	
+	// LMM related parameters
+	double l_min;
+	double l_max;
+	size_t n_region;
+	double pve_null;
+	double pheno_mean;
+	
+	// BSLMM MCMC related parameters
+	double h_min, h_max, h_scale;			//priors for h
+	double rho_min, rho_max, rho_scale;		//priors for rho
+	double logp_min, logp_max, logp_scale;		//priors for log(pi)
+	size_t s_min, s_max;			//minimum and maximum number of gammas
+	size_t w_step;					//number of warm up/burn in iterations
+	size_t s_step;					//number of sampling iterations
+	size_t r_pace;					//record pace
+	size_t w_pace;					//write pace
+	size_t n_accept;				//number of acceptance
+	size_t n_mh;					//number of MH steps within each iteration
+	double geo_mean;				//mean of the geometric distribution
+	long int randseed;
+	double trace_G;	
+	
+	HYPBSLMM cHyp_initial;
+
+	// 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
+	double time_UtZ;
+	double time_Omega;		//time spent on optimization iterations
+	double time_Proposal;        //time spent on constructing the proposal distribution for gamma (i.e. lmm or lm analysis)
+	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
+	
+	vector<SNPINFO> snpInfo;		//record SNP information
+	
+	// Not included in PARAM
+	gsl_rng *gsl_r; 
+	gsl_ran_discrete_t *gsl_t;	
+	map<size_t, size_t> mapRank2pos;	
+	
+	// Main Functions
+	void CopyFromParam (PARAM &cPar);
+	void CopyToParam (PARAM &cPar);
+	
+	void RidgeR(const gsl_matrix *U, const gsl_matrix *UtX, const gsl_vector *Uty, const gsl_vector *eval, const double lambda);
+	
+	void MCMC (const gsl_matrix *U, const gsl_matrix *UtX, const gsl_vector *Uty, const gsl_vector *K_eval, const gsl_vector *y);
+	void WriteLog ();
+	void WriteLR ();
+	void WriteBV (const gsl_vector *bv);
+	void WriteParam (vector<pair<double, double> > &beta_g, const gsl_vector *alpha, const size_t w);
+	void WriteParam (const gsl_vector *alpha);
+	void WriteResult (const int flag, const gsl_matrix *Result_hyp, const gsl_matrix *Result_gamma, const size_t w_col);
+	
+	//Subfunctions inside MCMC
+	void CalcPgamma (double *p_gammar);
+	
+	double CalcPveLM (const gsl_matrix *UtXgamma, const gsl_vector *Uty, const double sigma_a2);
+	void InitialMCMC (const gsl_matrix *UtX, const gsl_vector *Uty, vector<size_t> &rank_old, class HYPBSLMM &cHyp, vector<pair<size_t, double> > &pos_loglr);
+	double CalcPosterior (const gsl_vector *Uty, const gsl_vector *K_eval, gsl_vector *Utu, gsl_vector *alpha_prime, class HYPBSLMM &cHyp);
+	double CalcPosterior (const gsl_matrix *UtXgamma, const gsl_vector *Uty, const gsl_vector *K_eval, gsl_vector *UtXb, gsl_vector *Utu, gsl_vector *alpha_prime, gsl_vector *beta, class HYPBSLMM &cHyp);
+	void CalcCC_PVEnZ (const gsl_matrix *U, const gsl_vector *Utu, gsl_vector *z_hat, class HYPBSLMM &cHyp);
+	void CalcCC_PVEnZ (const gsl_matrix *U, const gsl_vector *UtXb, const gsl_vector *Utu, gsl_vector *z_hat, class HYPBSLMM &cHyp);
+	double CalcREMLE (const gsl_matrix *Utw, const gsl_vector *Uty, const gsl_vector *K_eval);
+	double CalcLR (const gsl_matrix *U, const gsl_matrix *UtX, const gsl_vector *Uty, const gsl_vector *K_eval, vector<pair<size_t, double> > &loglr_sort);		//calculate the maximum marginal likelihood ratio for each analyzed SNPs with gemma, use it to rank SNPs
+	void SampleZ (const gsl_vector *y, const gsl_vector *z_hat, gsl_vector *z);
+	double ProposeHnRho (const class HYPBSLMM &cHyp_old, class HYPBSLMM &cHyp_new, const size_t &repeat);
+	double ProposePi (const class HYPBSLMM &cHyp_old, class HYPBSLMM &cHyp_new, const size_t &repeat);
+	double ProposeGamma (const vector<size_t> &rank_old, vector<size_t> &rank_new, const double *p_gamma, const class HYPBSLMM &cHyp_old, class HYPBSLMM &cHyp_new, const size_t &repeat);
+	void SetXgamma (gsl_matrix *Xgamma, const gsl_matrix *X, vector<size_t> &rank);
+
+	void CalcXtX (const gsl_matrix *X_new, const gsl_vector *y, const size_t s_size, gsl_matrix *XtX_new, gsl_vector *Xty_new);
+	void SetXgamma (const gsl_matrix *X, const gsl_matrix *X_old, const gsl_matrix *XtX_old, const gsl_vector *Xty_old, const gsl_vector *y, const vector<size_t> &rank_old, const vector<size_t> &rank_new, gsl_matrix *X_new, gsl_matrix *XtX_new, gsl_vector *Xty_new);
+	double CalcPosterior (const double yty, class HYPBSLMM &cHyp);
+	double CalcPosterior (const gsl_matrix *Xgamma, const gsl_matrix *XtX, const gsl_vector *Xty, const double yty, const size_t s_size, gsl_vector *Xb, gsl_vector *beta, class HYPBSLMM &cHyp);
+	void CalcCC_PVEnZ (gsl_vector *z_hat, class HYPBSLMM &cHyp);
+	void CalcCC_PVEnZ (const gsl_vector *Xb, gsl_vector *z_hat, class HYPBSLMM &cHyp);
+	void MCMC (const gsl_matrix *X, const gsl_vector *y);
+	
+	//utility functions
+//	double vec_sum (gsl_vector *v);
+//	void vec_center (gsl_vector *v);
+//	double calc_var (gsl_vector *v);
+//	void calc_sigma (MCMC &cMcmc);
+//	bool comp_lr (pair<size_t, double> a, pair<size_t, double> b);
+};
+
+
+
+#endif
+
+