Output is on par with earlier gemma - only allele values are missing

4 files changed, 15 insertions, 13 deletions

diff --git a/src/lmm.cpp b/src/lmm.cpp
index 08b4d41..10e3e53 100644
--- a/src/lmm.cpp
+++ b/src/lmm.cpp
@@ -153,7 +153,7 @@ void LMM::WriteFiles() {
     outfile << scientific << setprecision(6);
 
     if (a_mode != M_LMM2) {
-      outfile << st.beta << "\t";
+      outfile << scientific << st.beta << "\t";
       outfile << st.se << "\t";
     }
 
@@ -2103,14 +2103,13 @@ void LMM::mdb_analyze(std::function< SnpNameValues2(size_t) >& fetch_snp,
     auto m = st.markerinfo;
     auto name = m.name;
     auto chr  = m.chr;
-    auto pos  = m.pos;
 
     outfile << chr << "\t";
     outfile << name << "\t";
-    outfile << pos << "\t";
-    auto n_miss = 0;
-    outfile << n_miss << "\t-\t-\t"; // n_miss column
-    // outfile << st.maf << "\t" ;
+    outfile << m.pos << "\t";
+    outfile << m.n_miss << "\t-\t-\t"; // n_miss column + allele1 + allele0
+    outfile << fixed << setprecision(3) << m.maf << "\t";
+    outfile << scientific << setprecision(6);
     if (a_mode != M_LMM2) {
       outfile << st.beta << "\t";
       outfile << st.se << "\t";
@@ -2253,9 +2252,9 @@ void LMM::mdb_calc_gwa(const gsl_matrix *U, const gsl_vector *eval,
 
       // compute maf and n_miss (NAs)
       size_t n_miss = 0; // count NAs: FIXME
-      double maf = compute_maf(ni_total, ni_test, n_miss, gs.data());
+      double maf = compute_maf(ni_total, ni_test, n_miss, gs.data(), indicator_idv);
 
-      markerinfo = MarkerInfo { .name=marker,.chr=chr,.pos=pos,.line_no=num };
+      markerinfo = MarkerInfo { .name=marker,.chr=chr,.pos=pos,.line_no=num,.n_miss=n_miss,.maf=maf };
 
       // cout << "!!!!" << size << marker << ": af" << maf << " " << gs[0] << "," << gs[1] << "," << gs[2] << "," << gs[3] << endl;
     }
diff --git a/src/lmm.h b/src/lmm.h
index 7643984..af79f5f 100644
--- a/src/lmm.h
+++ b/src/lmm.h
@@ -48,7 +48,8 @@ public:
 // typedef tuple< string, uint16_t, uint32_t, uint32_t > MarkerInfo; // name, chr, pos, line
 struct MarkerInfo {
   string name;
-  size_t chr, pos, line_no;
+  size_t chr, pos, line_no, n_miss;
+  double maf;
 } ;
 
 typedef vector<MarkerInfo> Markers;
diff --git a/src/mathfunc.cpp b/src/mathfunc.cpp
index f132e6e..c82872e 100644
--- a/src/mathfunc.cpp
+++ b/src/mathfunc.cpp
@@ -666,12 +666,14 @@ std::tuple<size_t, size_t, size_t> compute_ratio(size_t ni_total, const gsl_vect
   return {n_0,n_1,n_2};
 }
 
-double compute_maf(size_t ni_total, size_t ni_test, size_t n_miss, const double *gs) {
+double compute_maf(size_t ni_total, size_t ni_test, size_t n_miss, const double *gs, const vector<int> &indicator) {
   double maf = 0.0;
 
   for (size_t i = 0; i < ni_total; ++i) { // read genotypes
-    double geno = gs[i];
-    maf += geno; // 0..2 range
+    if (indicator[i]) {
+      double geno = gs[i];
+      maf += geno; // 0..2 range
+    }
   }
   maf /= 2.0 * (double)(ni_test - n_miss); // Assumption is that geno value is between 0 and 2. FIXME
   return maf;
diff --git a/src/mathfunc.h b/src/mathfunc.h
index b408ea7..f7c9e6e 100644
--- a/src/mathfunc.h
+++ b/src/mathfunc.h
@@ -79,7 +79,7 @@ unsigned char Double02ToUchar(const double dosage);
 //                          const size_t i_row, Eigen::VectorXd &x_row);
 
 std::tuple<size_t, size_t, size_t> compute_ratio(size_t ni_total, const gsl_vector *gs);
-double compute_maf(size_t ni_total, size_t ni_test, size_t n_miss, const double *gs);
+double compute_maf(size_t ni_total, size_t ni_test, size_t n_miss, const double *gs, const vector<int> &indicator);
 
 
 #endif