about summary refs log tree commit diff
diff options
context:
space:
mode:
authorPjotr Prins2025-12-09 09:18:48 +0100
committerPjotr Prins2025-12-09 09:18:48 +0100
commit6a18d7ffd730dea629ba5198ecd13b00e946a175 (patch)
treea0589af2f6be9f13e155d0ca69fc6519d5888252
parent64ade49ee212fca39087509505bd95d07df46c11 (diff)
downloadpangemma-6a18d7ffd730dea629ba5198ecd13b00e946a175.tar.gz
Add header to output
-rw-r--r--src/lmm.cpp47
1 files changed, 38 insertions, 9 deletions
diff --git a/src/lmm.cpp b/src/lmm.cpp
index 2b730f3..9fe6bbb 100644
--- a/src/lmm.cpp
+++ b/src/lmm.cpp
@@ -1981,7 +1981,6 @@ void LMM::mdb_analyze(std::function< SnpNameValues2(size_t) >& fetch_snp,
         CalcLambda('L', param1, l_min, l_max, n_region, lambda_mle, logl_H1);
         p_lrt = gsl_cdf_chisq_Q(2.0 * (logl_H1 - logl_mle_H0), 1);
       }
-
       time_opt += (clock() - time_start) / (double(CLOCKS_PER_SEC) * 60.0);
 
       auto markerinfo = markers[i];
@@ -1991,14 +1990,6 @@ void LMM::mdb_analyze(std::function< SnpNameValues2(size_t) >& fetch_snp,
     }
   };
 
-  /*
-  const auto num_markers = indicator_snp.size();
-  enforce_msg(num_markers > 0,"Zero SNPs to process - data corrupt?");
-  if (num_markers < 50) {
-    cerr << num_markers << " SNPs" << endl;
-    warning_msg("very few SNPs processed");
-  }
-  */
   const size_t progress_step = (num_markers/50>d_pace ? num_markers/50 : d_pace);
   Markers markers;
 
@@ -2087,6 +2078,44 @@ void LMM::mdb_analyze(std::function< SnpNameValues2(size_t) >& fetch_snp,
     return;
   }
 
+  outfile << "chr" << "\t"
+          << "rs" << "\t"
+          << "ps" << "\t"
+          << "n_miss" << "\t"
+          << "allele1" << "\t"
+          << "allele0" << "\t"
+          << "af" << "\t";
+
+  if (a_mode != M_LMM2) {
+    outfile << "beta" << "\t";
+    outfile << "se" << "\t";
+  }
+
+  if (a_mode != M_LMM3 && a_mode != M_LMM9)
+    outfile << "logl_H1" << "\t";
+
+  switch(a_mode) {
+  case M_LMM1:
+    outfile << "l_remle" << "\t"
+            << "p_wald" << endl;
+    break;
+  case M_LMM2:
+  case M_LMM9:
+    outfile << "l_mle" << "\t"
+            << "p_lrt" << endl;
+    break;
+  case M_LMM3:
+    outfile << "p_score" << endl;
+    break;
+  case M_LMM4:
+    outfile << "l_remle" << "\t"
+            << "l_mle" << "\t"
+            << "p_wald" << "\t"
+            << "p_lrt" << "\t"
+            << "p_score" << endl;
+    break;
+  }
+
   auto sumstats = [&] (SUMSTAT2 st) {
     outfile << scientific << setprecision(6);
     auto m = st.markerinfo;