review update

Author: vasudeva8

htslib/vcf.h

@@ -1523,15 +1523,11 @@ int bcf_format_gt_v2(const bcf_hdr_t *hdr, bcf_fmt_t *fmt, int isample,
  *  update44phasing - converts GT to/from v4.4 way representation
  *  @param h - bcf header, to get version
  *  @param v - pointer to bcf data
- *  @param setreset - whether to set or reset
- *  Returns 0 on success and -1 on failure
- *  For data read, to be converted to v44, setreset to be 1. For data write, to
- *  be converted to v < v44, setreset to be 0.
+ *  Returns 0 on success and 1 on failure
  *  If the version in header is >= 4.4, no change is made. Otherwise 1st phasing
  *  is set if there are no other unphased ones.
  */
-HTSLIB_EXPORT int update44phasing(bcf_hdr_t *h, bcf1_t *b, int setreset)
- HTS_RESULT_USED;
+HTSLIB_EXPORT int update44phasing(bcf_hdr_t *h, bcf1_t *b) HTS_RESULT_USED;
 
 static inline int bcf_format_gt(bcf_fmt_t *fmt, int isample, kstring_t *str)
 {

synced_bcf_reader.c

@@ -709,7 +709,7 @@ static int _reader_fill_buffer(bcf_srs_t *files, bcf_sr_t *reader)
             if ( ret < -1 ) files->errnum = bcf_read_error;
             if ( ret < 0 ) break; // no more lines or an error
             //set phasing of 1st value as in vcf v44
-            if ((ret = update44phasing(reader->header, reader->buffer[reader->nbuffer+1], 1))) {
+            if ((ret = update44phasing(reader->header, reader->buffer[reader->nbuffer+1]))) {
                 files->errnum = bcf_read_error;
                 break;
             }

vcf.c

@@ -1819,64 +1819,49 @@ static void bcf_record_check_err(const bcf_hdr_t *hdr, bcf1_t *rec,
  *  @param p - pointer to phase value array
  *  @param end - end of array
  *  @param q - pointer to consumed data
- *  @param set - whether to set (while read) or reset (for write)
  *  @param samples - no. of samples in array
  *  @param ploidy - no. of phasing values per sample
  *  @param type - value type (one of BCF_BT_...)
  *  Returns 0 on success and 1 on failure
  */
-static int updatephasing(uint8_t *p, uint8_t *end, uint8_t **q, int set, int samples, int ploidy, int type)
+static int updatephasing(uint8_t *p, uint8_t *end, uint8_t **q, int samples, int ploidy, int type)
 {
-    int j, k, upd = set ? 1 : -1;
+    int j, k;
     size_t bytes;
-    for (j = 0; j < samples; ++j) {   //for each sample
-        int anyunphased = 0;
-        uint8_t *ptr1 = p;
-        int32_t al1 = 0, val;
-        for (k = 0; k < ploidy; ++k) { //in each ploidy
-            switch(type) {
-                case BCF_BT_INT8:
-                    val = *(uint8_t*)p;
-                break;
-                case BCF_BT_INT16:
-                    val = *(uint16_t*)p;
-                break;
-                case BCF_BT_INT32:
-                    val = *(uint32_t*)p;
-                break;
-                //wont have anything bigger than 32bit for GT
-                default: //invalid
-                    return 1;
-                break;
-            }
-            if (!k) al1 = val;
-            else if (!(val & 1)) {
-                anyunphased = 1;
-            }
-            //get to next phasing or skip the rest for this sample
-            bytes = (anyunphased ? ploidy - k : 1) << bcf_type_shift[type];
-            if (end - p < bytes)
-                return 1;
-            p += bytes;
-            if (anyunphased) {
-                break;  //no further check required
-            }
-        }
-        if (!anyunphased && al1 > 1) {  //no other unphased
+
+    #define NOSWITCHFOR(type_t, rawtype_t,convert, reconvert, p, end, samples, ploidy, consumed) \
+     for (j = 0; j < samples; ++j) {   /*for each sample */ \
+        int anyunphased = 0; \
+        uint8_t *ptr1 = p; \
+        rawtype_t al1 = 0, val; \
+        for (k = 0; k < ploidy; ++k) { /*in each ploidy */ \
+            val = convert(p); \
+            if (!k) al1 = val; \
+            else if (!(val & 1)) anyunphased = 1; \
+            /*get to next phasing or skip the rest for this sample*/ \
+            *consumed = ((anyunphased || al1 & 1) ? ploidy - k : 1) << bcf_type_shift[type_t]; \
+            if (end - p < *consumed) return 1; \
+            p += *consumed; \
+            if (anyunphased || al1 & 1) break; \
+        } \
+        if (!anyunphased && al1 > 1) {  /*no other unphased*/ \
             /*set phased on read or make unphased on write as upto 4.3 1st
-            phasing is not described explicitly and has to be inferred*/
-            switch(type) {
-                case BCF_BT_INT8:
-                    *(uint8_t*)ptr1 += upd;
-                break;
-                case BCF_BT_INT16:
-                    *(uint16_t*)ptr1 += upd;
-                break;
-                case BCF_BT_INT32:
-                    *(uint32_t*)ptr1 += upd;
-                break;
-            }
-        }
+            phasing is not described explicitly and has to be inferred*/ \
+            al1 |= 1; \
+            reconvert(convert((uint8_t*)&al1), ptr1); \
+        } \
+    }
+    #define i8_to_le(v, buf) *((uint8_t*)buf) = v;
+    switch(type) {
+        case BCF_BT_INT8:
+        NOSWITCHFOR(BCF_BT_INT8, uint8_t, le_to_i8, i8_to_le, p, end, samples, ploidy, &bytes);
+        break;
+        case BCF_BT_INT16:
+        NOSWITCHFOR(BCF_BT_INT16, uint16_t, le_to_i16, i16_to_le, p, end, samples, ploidy, &bytes);
+        break;
+        case BCF_BT_INT32:
+        NOSWITCHFOR(BCF_BT_INT32, uint32_t, le_to_i32, i32_to_le, p, end, samples, ploidy, &bytes);
+        break;
     }
     *q = p;
     return 0;
@@ -1886,14 +1871,11 @@ static int updatephasing(uint8_t *p, uint8_t *end, uint8_t **q, int set, int sam
  *  update44phasing - converts GT to/from v4.4 way representation
  *  @param h - bcf header, to get version
  *  @param v - pointer to bcf data
- *  @param setreset - whether to set or reset
- *  Returns 0 on success and -1 on failure
- *  For data read, to be converted to v44, setreset to be 1. For data write, to
- *  be converted to v < v44, setreset to be 0.
+ *  Returns 0 on success and 1 on failure
  *  If the version in header is >= 4.4, no change is made. Otherwise 1st phasing
  *  is set if there are no other unphased ones.
  */
-HTSLIB_EXPORT int update44phasing(bcf_hdr_t *h, bcf1_t *b, int setreset)
+HTSLIB_EXPORT int update44phasing(bcf_hdr_t *h, bcf1_t *b)
 {
     int i, idgt = -1, ver = VCF_DEF, num, type;
     uint8_t *ptr = NULL, *end = NULL;
@@ -1938,7 +1920,7 @@ HTSLIB_EXPORT int update44phasing(bcf_hdr_t *h, bcf1_t *b, int setreset)
     }
     if (ptr) {
         //with GT and v < v44, need phase conversion
-        if (updatephasing(ptr, end, &ptr, setreset, b->n_sample, num, type)) return 1;
+        if (updatephasing(ptr, end, &ptr, b->n_sample, num, type)) return 1;
     }
     return 0;
 }
@@ -2076,8 +2058,8 @@ static int bcf_record_check(const bcf_hdr_t *hdr, bcf1_t *rec) {
         }
         if (ver < VCF44 && idgt != -1 && idgt == key) {
             //with GT and v < v44, need phase conversion
-            if (updatephasing(ptr, end, &ptr, 1, rec->n_sample, num, type)) {
-                err |= BCF_ERR_TAG_INVALID;
+            if (updatephasing(ptr, end, &ptr, rec->n_sample, num, type)) {
+               err |= BCF_ERR_TAG_INVALID;
             }
         } else {
             bytes = ((size_t) num << bcf_type_shift[type]) * rec->n_sample;
@@ -2442,11 +2424,6 @@ int bcf_write(htsFile *hfp, bcf_hdr_t *h, bcf1_t *v)
         return -1;
     }
 
-    if (update44phasing(h, v, 0)) { //reset phasing update made after read
-        hts_log_error("Failed to set prorper phasing at %s:%"PRIhts_pos"", bcf_seqname_safe(h,v), v->pos+1);
-        return -1;
-    }
-
     BGZF *fp = hfp->fp.bgzf;
     uint8_t x[32];
     u32_to_le(v->shared.l + 24, x); // to include six 32-bit integers