debian/changelog

+plink1.9 (1.90~b6.5-180913-1) unstable; urgency=medium
+
+  * New upstream release.
+  * Bump Standards-Version: 4.2.1 (no changes needed).
+
+ -- Dylan Aïssi <bob.dybian@gmail.com>  Sun, 16 Sep 2018 15:45:33 +0200
+
 plink1.9 (1.90~b6.4-180807-1) unstable; urgency=medium
 
   * New upstream release.

debian/control

@@ -8,7 +8,7 @@ Build-Depends: debhelper (>= 11~),
                libatlas-base-dev,
                liblapack-dev,
                zlib1g-dev
-Standards-Version: 4.2.0
+Standards-Version: 4.2.1
 Vcs-Browser: https://salsa.debian.org/med-team/plink1-9
 Vcs-Git: https://salsa.debian.org/med-team/plink1-9.git
 Homepage: https://www.cog-genomics.org/plink/1.9/

debian/upstream.docs/upstream.changelog

 # Copy/Paste from https://www.cog-genomics.org/plink/1.9/
 
-7 Aug 2018: --assoc set-test bugfix.
+13 Sep 2018: Improved logistic regression convergence-failure detection (affects case/control --epistasis as well as --logistic). Fixed --R memory leak.
+
+7 Aug: --assoc set-test bugfix.
 
 17 Jul: Fixed uninitialized --make-set-border value bug.
 

plink.c

@@ -93,7 +93,7 @@
 static const char ver_str[] =
 #ifdef STABLE_BUILD
-  "PLINK v1.90b6.4"
+  "PLINK v1.90b6.5"
 #else
   "PLINK v1.90p"
 #endif
@@ -105,7 +105,7 @@ static const char ver_str[] =
 #else
   " 32-bit"
 #endif
-  " (7 Aug 2018)";
+  " (13 Sep 2018)";
 static const char ver_str2[] =
   // include leading space if day < 10, so character length stays the same
   ""

plink_calc.c

@@ -3028,7 +3028,7 @@ int32_t groupdist_calc(pthread_t* threads, uint32_t unfiltered_sample_ct, uintpt
     dxx = g_reg_tot_x / dww;
     dhh_sd = sqrt((dhh_ssq / dww - dxx * dxx) / (dww - 1.0));
   }
-  if (!(g_case_ct * g_ctrl_ct)) {
+  if (!(g_case_ct && g_ctrl_ct)) {
     dyy = 0.0;
     dhl_sd = 0.0;
   } else {

plink_data.c

@@ -6184,7 +6184,7 @@ int32_t ped_to_bed(char* pedname, char* mapname, char* outname, char* outname_en
       if (bigstack_left() >= marker_ct * sample_ct4) {
 	markers_per_pass = marker_ct;
 	sprintf(g_logbuf, "Performing single-pass .bed write (%" PRIuPTR " variant%s, %" PRIuPTR " %s).\n", marker_ct, (marker_ct == 1)? "" : "s", sample_ct, species_str(sample_ct));
-	pass_ct = (marker_ct * sample_ct4)? 1 : 0;
+	pass_ct = (marker_ct && sample_ct4)? 1 : 0;
       } else {
 	if (!map_is_unsorted) {
 	  if (bigstack_alloc_ll(sample_ct, &line_starts)) {

plink_glm.c

@@ -1658,12 +1658,29 @@ uint32_t logistic_regression(uint32_t sample_ct, uint32_t param_ct, float* vv, f
 	return 1;
       }
       if (iteration >= 15) {
+        // quasi-bugfix (11 Sep 2018): if fabsf(any coefficient) > 8e3, this is
+        // almost certainly a form of convergence failure that didn't get
+        // caught by the (delta_coef > 20.0) check due to a precision quirk.
+        // (8e3 threshold ~= 1e-4 * 2^23, since floats have 23 bits of
+        // precision)
+        for (param_idx = 0; param_idx < param_ct; param_idx++) {
+          if (fabsf(coef[param_idx]) > 8e3) {
+            return 1;
+          }
+        }
 	return 0;
       }
     }
     // Pons reported that 1.1e-3 was dangerous, so I agree with the decision to
     // tighten this threshold from 1e-3 to 1e-4.
     if (delta_coef < 1e-4) {
+      // Be more conservative in throwing out results when we don't hit the
+      // iteration limit.
+      for (param_idx = 0; param_idx < param_ct; param_idx++) {
+        if (fabsf(coef[param_idx]) > 6e4) {
+          return 1;
+        }
+      }
       return 0;
     }
   }

plink_ld.c

@@ -3233,6 +3233,7 @@ void fepi_counts_to_joint_effects_stats(uint32_t group_ct, uint32_t* counts, dou
   uint32_t uii;
   uint32_t ujj;
   uint32_t ukk;
+  tot_inv_v[0] = 0.0;  // gcc7 maybe-uninitialized warning
   dptr = dcounts;
   if (counts[0] && counts[1] && counts[2] && counts[3] && counts[4] && counts[5] && counts[6] && counts[7] && counts[8] && ((group_ct == 1) || (counts[9] && counts[10] && counts[11] && counts[12] && counts[13] && counts[14] && counts[15] && counts[16] && counts[17]))) {
     for (uii = 0; uii < group_ct; uii++) {
@@ -4712,11 +4713,15 @@ THREAD_RET_TYPE epi_logistic_thread(void* arg) {
 	  cur_word2 = cur_word2 & (cur_word2 >> 1);
 	  cur_sample_ct -= popcount2_long((cur_word1 | cur_word2) & FIVEMASK);
 	}
+        unsigned char geno_pair_present[12];
 	if (cur_sample_ct <= 4) {
 	  goto epi_logistic_thread_regression_fail;
 	}
 	// 2. now populate covariate-major matrix with 16-byte-aligned,
 	//    trailing-entries-zeroed rows
+        // quasi-bugfix (13 Sep 2018): reliably detect when this matrix is not
+        // of full rank, and skip the regression in that case.
+        memset(geno_pair_present, 0, 12);
 	cur_sample_cta4 = round_up_pow2(cur_sample_ct, 4);
 	for (widx = 0; widx < pheno_nm_ctl2; widx++) {
 	  sample_idx = widx * BITCT2;
@@ -4731,6 +4736,7 @@ THREAD_RET_TYPE epi_logistic_thread(void* arg) {
             uljj = cur_word2 & (3 * ONELU);
 	    if ((ulii != 3) && (uljj != 3)) {
               *fptr = 1.0;
+              geno_pair_present[ulii + uljj * 4] = 1;
 	      fxx = (float)((intptr_t)ulii);
 	      fyy = (float)((intptr_t)uljj);
 	      // maybe this is faster with continuous writes instead of
@@ -4745,6 +4751,16 @@ THREAD_RET_TYPE epi_logistic_thread(void* arg) {
             cur_word2 >>= 2;
 	  }
 	}
+        if (!geno_pair_present[5]) {
+          // not full rank if any 2x2 square in the 3x3 contingency table is
+          // empty.
+          if (((!geno_pair_present[0]) && (!geno_pair_present[1]) && (!geno_pair_present[4])) ||
+              ((!geno_pair_present[1]) && (!geno_pair_present[2]) && (!geno_pair_present[6])) ||
+              ((!geno_pair_present[4]) && (!geno_pair_present[8]) && (!geno_pair_present[9])) ||
+              ((!geno_pair_present[6]) && (!geno_pair_present[9]) && (!geno_pair_present[10]))) {
+            goto epi_logistic_thread_regression_fail;
+          }
+        }
 	if (cur_sample_ct < cur_sample_cta4) {
 	  loop_end = cur_sample_cta4 - cur_sample_ct;
 	  fill_float_zero(loop_end, fptr);

plink_rserve.c

@@ -63,6 +63,7 @@ int32_t rserve_call(char* rplugin_fname, char* rplugin_host_or_socket, int32_t r
   Rinteger* r_l;
   Rinteger* r_g;
   Rdouble* r_data;
+  Rexp* r_eval_dummy;  // bugfix (11 Sep 2018): must delete rc->eval return-val
   double* pheno_d2;
   double* covar_d2 = nullptr;
   double* dptr;
@@ -81,6 +82,12 @@ int32_t rserve_call(char* rplugin_fname, char* rplugin_host_or_socket, int32_t r
   uint32_t cur_data_len;
   uint32_t uii;
   int32_t ii;
+  if (!pheno_nm_ct) {
+    // bugfix (11 Sep 2018): this was segfaulting instead of printing an error
+    // message
+    logerrprint("Error: --R only processes samples with nonmissing phenotype values, and all\nphenotype values are missing.\n");
+    goto rserve_call_ret_1;
+  }
   if (fopen_checked(rplugin_fname, "r", &infile)) {
     goto rserve_call_ret_OPEN_FAIL;
   }
@@ -198,13 +205,19 @@ int32_t rserve_call(char* rplugin_fname, char* rplugin_host_or_socket, int32_t r
     rc->assign("n", r_n);
     rc->assign("PHENO", r_p);
     rc->assign("CLUSTER", r_s);
-    rc->eval("CLUSTER[CLUSTER==-1] <- NA");
+    r_eval_dummy = rc->eval("CLUSTER[CLUSTER==-1] <- NA");
+    if (r_eval_dummy) {
+      delete r_eval_dummy;
+    }
     if (covar_ct) {
       r_cov = new Rdouble(covar_d2, pheno_nm_ct * covar_ct);
       rc->assign("c", r_cov);
-      rc->eval("COVAR<-matrix(c,nrow=n,byrow=T)");
+      r_eval_dummy = rc->eval("COVAR<-matrix(c,nrow=n,byrow=T)");
     } else {
-      rc->eval("COVAR<-NA");
+      r_eval_dummy = rc->eval("COVAR<-NA");
+    }
+    if (r_eval_dummy) {
+      delete r_eval_dummy;
     }
     memcpy(outname_end, ".auto.R", 8);
   } else {
@@ -314,11 +327,20 @@ int32_t rserve_call(char* rplugin_fname, char* rplugin_host_or_socket, int32_t r
       r_g = new Rinteger(geno_int_buf, block_size * pheno_nm_ct);
       rc->assign("l", r_l);
       rc->assign("g", r_g);
-      rc->eval("GENO<-matrix(g,nrow=n,byrow=T)");
-      rc->eval("GENO[GENO==-1] <- NA");
+      r_eval_dummy = rc->eval("GENO<-matrix(g,nrow=n,byrow=T)");
+      if (r_eval_dummy) {
+        delete r_eval_dummy;
+      }
+      r_eval_dummy = rc->eval("GENO[GENO==-1] <- NA");
+      if (r_eval_dummy) {
+        delete r_eval_dummy;
+      }
       delete r_l;
       delete r_g;
-      rc->eval(inbuf_start);
+      r_eval_dummy = rc->eval(inbuf_start);
+      if (r_eval_dummy) {
+        delete r_eval_dummy;
+      }
       r_data = (Rdouble*)rc->eval("Rplink(PHENO,GENO,CLUSTER,COVAR)");
       if (r_data) {
 	dptr = r_data->doubleArray();
@@ -463,6 +485,7 @@ int32_t rserve_call(char* rplugin_fname, char* rplugin_host_or_socket, int32_t r
     retval = RET_NETWORK;
     break;
   }
+ rserve_call_ret_1:
   fclose_cond(infile);
   fclose_cond(outfile);
   bigstack_reset(bigstack_mark);