Merge pull request #9 from MRCIEU/openmp

Openmp

sim/sim4.R

@@ -19,6 +19,18 @@ for (t in c(1,2,4,8)){
         # simulate covariates
         data <- data.frame(
             S = paste0("S", seq(1, n_obs)),
+            Age = rnorm(n_obs),
+            Sex = rbinom(n_obs, 2, .5),
+            PC1 = rnorm(n_obs),
+            PC2 = rnorm(n_obs),
+            PC3 = rnorm(n_obs),
+            PC4 = rnorm(n_obs),
+            PC5 = rnorm(n_obs),
+            PC6 = rnorm(n_obs),
+            PC7 = rnorm(n_obs),
+            PC8 = rnorm(n_obs),
+            PC9 = rnorm(n_obs),
+            PC10 = rnorm(n_obs),
             Y = rnorm(n_obs),
             stringsAsFactors=F
         )

src/Model.cpp

@@ -29,6 +29,7 @@ void Model::run() {
   std::vector<std::string> alleles;
   std::vector<std::vector<double>> probs;
   std::vector<double> dosages;
+  std::string a0, a1;
   spdlog::info("Using {} thread(s)", _threads);
   omp_set_num_threads(_threads);
 
@@ -62,7 +63,7 @@ void Model::run() {
     file << "chr\tpos\trsid\toa\tea\tn\teaf\tbeta\tse\tt\tp\ttheta\tphi_x1\tse_x1\tphi_x2\tse_x2\tphi_f\tphi_p"
          << std::endl;
     file.flush();
-#pragma omp parallel default(none) shared(file, X1, X2, y) private(chromosome, position, rsid, alleles, probs, dosages)
+#pragma omp parallel default(none) shared(file, X1, X2, y) private(chromosome, position, rsid, alleles, probs, dosages, a0, a1)
     {
 #pragma omp master
       // Read variant-by-variant
@@ -74,6 +75,14 @@ void Model::run() {
           continue;
         }
 
+        if (_flip) {
+          a0 = alleles[1];
+          a1 = alleles[0];
+        } else {
+          a0 = alleles[0];
+          a1 = alleles[1];
+        }
+
         // convert probabilities to dosage values
         _bgen_parser.read_probs(&probs);
 
@@ -98,7 +107,11 @@ void Model::run() {
           if ((prob[0] == -1 && prob[1] == -1 && prob[2] == -1) || (prob[0] == 0 && prob[1] == 0 && prob[2] == 0)) {
             dosages.push_back(-1);
           } else {
-            dosages.push_back(prob[1] + (2 * prob[0]));
+            if (_flip) {
+              dosages.push_back(prob[1] + (2 * prob[2]));
+            } else {
+              dosages.push_back(prob[1] + (2 * prob[0]));
+            }
           }
         }
 
@@ -115,8 +128,8 @@ void Model::run() {
               res = fit(chromosome,
                         position,
                         rsid,
-                        alleles[0],
-                        alleles[1],
+                        a0,
+                        a1,
                         dosages,
                         _non_null_idx,
                         X1,
@@ -219,14 +232,29 @@ void Model::delta_method(const Eigen::VectorXd &ss_fit,
                          double &s2_dummy) {
   spdlog::debug("Preparing partial derivatives for deltamethod, thread = {}", omp_get_thread_num());
   const double pi = boost::math::constants::pi<double>();
-  Eigen::MatrixXd grad1 = Eigen::MatrixXd(3, 1);
-  grad1(0, 0) = 0;
-  grad1(1, 0) = (2 * ss_fit(0, 0) + 2 * ss_fit(1, 0)) / (2 / pi);
-  grad1(2, 0) = 0;
-  Eigen::MatrixXd grad2 = Eigen::MatrixXd(3, 1);
-  grad2(0, 0) = 0;
-  grad2(1, 0) = 0;
-  grad2(2, 0) = (2 * ss_fit(0, 0) + 2 * ss_fit(2, 0)) / (2 / pi);
+
+  Eigen::MatrixXd grad1 = Eigen::MatrixXd(ss_fit.rows(), 1);
+  grad1(0, 0) = 0; // intercept
+  grad1(1, 0) = (2 * ss_fit(0, 0) + 2 * ss_fit(1, 0)) / (2 / pi); // SNP=1
+  grad1(2, 0) = 0; // SNP=2
+
+  Eigen::MatrixXd grad2 = Eigen::MatrixXd(ss_fit.rows(), 1);
+  grad2(0, 0) = 0; // intercept
+  grad2(1, 0) = 0; // SNP=1
+  grad2(2, 0) = (2 * ss_fit(0, 0) + 2 * ss_fit(2, 0)) / (2 / pi); // SNP=2
+
+  // fill in covariates as zeros
+  for (unsigned j = 3; j < ss_fit.rows(); j++) {
+    grad1(j, 0) = 0;
+    grad2(j, 0) = 0;
+  }
+
+  if (grad1.rows() != ss_fit.rows()) {
+    throw std::runtime_error("Number of observations for G1 and X differ");
+  }
+  if (grad2.rows() != ss_fit.rows()) {
+    throw std::runtime_error("Number of observations for G2 and X differ");
+  }
 
   spdlog::debug("Estimating HC SE_1 using the deltamethod, thread = {}", omp_get_thread_num());
   Eigen::MatrixXd s1_hc0 = grad1.transpose() * hc_vcov * grad1;

src/Model.h

@@ -28,13 +28,16 @@ class Model {
         std::set<unsigned> &non_null_idx,
         std::string &output_file,
         int threads,
-        double maf_threshold)
+        double maf_threshold,
+        bool flip)
       : _phenotype_file(phenotype_file),
         _bgen_parser(bgen_parser),
         _non_null_idx(non_null_idx),
         _output_file(output_file),
         _threads(threads),
-        _maf_threshold(maf_threshold) {}
+        _maf_threshold(maf_threshold),
+        _flip(flip)
+        {}
   void run();
   static Result fit(std::string &chromosome,
                     uint32_t position,
@@ -53,6 +56,7 @@ class Model {
   std::set<unsigned> &_non_null_idx;
   int _threads;
   double _maf_threshold;
+  bool _flip;
   std::string &_output_file;
   static void first_stage_ols(
       const Eigen::MatrixXd &X_complete1,

src/main.cpp

@@ -18,7 +18,7 @@ bool file_exists(const std::string &name) {
 }
 
 int main(int argc, char **argv) {
-  static std::string VERSION = "v1.2.2";
+  static std::string VERSION = "v1.2.3";
   static std::string PROGRAM_NAME = "varGWAS";
   spdlog::cfg::load_env_levels();
   static bool no_args = (argc == 1);
@@ -38,6 +38,7 @@ int main(int argc, char **argv) {
       ("p,phenotype", "Column name for phenotype", cxxopts::value<std::string>())
       ("i,id", "Column name for genotype identifier", cxxopts::value<std::string>())
       ("m,maf", "Filter out variants with a MAF below this threshold", cxxopts::value<double>())
+      ("f,flip", "Flip alleles")
       ("h,help", "Print usage")
       ("t,threads",
        "Number of threads",
@@ -110,7 +111,9 @@ int main(int argc, char **argv) {
   if (result.count("maf") == 1) {
     maf_threshold = result["maf"].as<double>();
   }
-  spdlog::debug("MAF threshold {}", maf_threshold);
+  spdlog::debug("MAF threshold: {}", maf_threshold);
+  bool flip = result.count("flip") == 1;
+  spdlog::debug("Flip: {}", flip);
 
   // check files exist
   if (!file_exists(variable_file)) {
@@ -141,7 +144,7 @@ int main(int argc, char **argv) {
     std::set<unsigned> non_null_idx = phenotype_file.join(samples);
 
     // Perform locus association tests & write to file
-    vargwas::Model model(phenotype_file, bgen_parser, non_null_idx, output_file, threads, maf_threshold);
+    vargwas::Model model(phenotype_file, bgen_parser, non_null_idx, output_file, threads, maf_threshold, flip);
     model.run();
 
     spdlog::info("Analysis complete!");