converting bcf 1st phase data as in v44

Author: vasudeva8

htslib/vcf.h

@@ -1518,6 +1518,21 @@ HTSLIB_EXPORT
 int bcf_format_gt_v2(const bcf_hdr_t *hdr, bcf_fmt_t *fmt, int isample,
     kstring_t *str) HTS_RESULT_USED;
 
+
+/**
+ *  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.
+ *  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;
+
 static inline int bcf_format_gt(bcf_fmt_t *fmt, int isample, kstring_t *str)
 {
     return bcf_format_gt_v2(NULL, fmt, isample, str);

synced_bcf_reader.c

@@ -708,6 +708,11 @@ static int _reader_fill_buffer(bcf_srs_t *files, bcf_sr_t *reader)
             ret = bcf_itr_next(reader->file, reader->itr, reader->buffer[reader->nbuffer+1]);
             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))) {
+                files->errnum = bcf_read_error;
+                break;
+            }
             bcf_subset_format(reader->header,reader->buffer[reader->nbuffer+1]);
         }
 

test/test.pl

@@ -369,9 +369,9 @@ sub failed
     if ( defined $reason ) { print STDERR "\t$reason\n"; }
     print STDERR ".. failed ...\n\n";
     STDERR->flush();
-    if ($$opts{fail_fast}) {
-      die "\n";
-    }
+    #if ($$opts{fail_fast}) {
+    #  die "\n";
+    #}
 }
 sub passed
 {

vcf.c

@@ -1814,6 +1814,135 @@ static void bcf_record_check_err(const bcf_hdr_t *hdr, bcf1_t *rec,
     (*reports)++;
 }
 
+/**
+ *  updatephasing - updates 1st phasing based on other phasing status
+ *  @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)
+{
+    int j, k, upd = set ? 1 : -1;
+    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
+            /*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;
+            }
+        }
+    }
+    *q = p;
+    return 0;
+}
+
+/**
+ *  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.
+ *  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)
+{
+    int i, idgt = -1, ver = VCF_DEF, num, type;
+    uint8_t *ptr = NULL, *end = NULL;
+    if (!b) return 1;
+
+    ver = bcf_get_version(h, "");
+    idgt = h ? bcf_hdr_id2int(h, BCF_DT_ID, "GT") : -1;
+    if (ver >= VCF44 || idgt == -1) return 0;   //no change required
+
+    if (b->unpacked & BCF_UN_FMT) { //unpacked, get from decoded data
+        for (i=0; i<b->n_fmt; i++)
+        {
+            if ( b->d.fmt[i].id == idgt ) {
+                ptr = b->d.fmt[i].p;
+                end = ptr + b->d.fmt[i].p_len;
+                num = b->d.fmt[i].n;
+                type = b->d.fmt[i].type;
+                break;
+            }
+        }
+    } else {    //get from indiv.s binary stream
+        ptr = (uint8_t *) b->indiv.s;
+        end = ptr + b->indiv.l;
+        int found = 0;
+        for (i = 0; i < b->n_fmt; ++i) {
+            int32_t key = -1;
+            if (bcf_dec_typed_int1_safe(ptr, end, &ptr, &key) != 0) return 1;
+            if (bcf_dec_size_safe(ptr, end, &ptr, &num, &type) != 0) return 1;
+            if (type > BCF_BT_CHAR) return 1;   //invalid type
+            if (idgt == key) {
+                found = 1;
+                break;
+            } else {    //skip and check next
+                size_t bytes = ((size_t) num << bcf_type_shift[type]) * b->n_sample;
+                if (end - ptr < bytes) return 1;
+                ptr += bytes;
+            }
+        }
+        if (!found) {
+            ptr = end = NULL;
+        }
+    }
+    if (ptr) {
+        //with GT and v < v44, need phase conversion
+        if (updatephasing(ptr, end, &ptr, setreset, b->n_sample, num, type)) return 1;
+    }
+    return 0;
+}
+
 static int bcf_record_check(const bcf_hdr_t *hdr, bcf1_t *rec) {
     uint8_t *ptr, *end;
     size_t bytes;
@@ -1832,6 +1961,10 @@ static int bcf_record_check(const bcf_hdr_t *hdr, bcf1_t *rec) {
                                     (1 << BCF_BT_FLOAT) |
                                     (1 << BCF_BT_CHAR));
     int32_t max_id = hdr ? hdr->n[BCF_DT_ID] : 0;
+    int idgt = -1, ver = VCF_DEF;
+
+    ver = bcf_get_version(hdr, NULL);
+    idgt = hdr ? bcf_hdr_id2int(hdr, BCF_DT_ID, "GT") : -1;
 
     // Check for valid contig ID
     if (rec->rid < 0
@@ -1941,9 +2074,16 @@ static int bcf_record_check(const bcf_hdr_t *hdr, bcf1_t *rec) {
             bcf_record_check_err(hdr, rec, "type", &reports, type);
             err |= BCF_ERR_TAG_INVALID;
         }
-        bytes = ((size_t) num << bcf_type_shift[type]) * rec->n_sample;
-        if (end - ptr < bytes) goto bad_indiv;
-        ptr += bytes;
+        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;
+            }
+        } else {
+            bytes = ((size_t) num << bcf_type_shift[type]) * rec->n_sample;
+            if (end - ptr < bytes) goto bad_indiv;
+            ptr += bytes;
+        }
     }
 
     if (!err && rec->rlen < 0) {
@@ -2032,6 +2172,8 @@ int bcf_readrec(BGZF *fp, void *null, void *vv, int *tid, hts_pos_t *beg, hts_po
     if (ret == 0) ret = bcf_record_check(NULL, v);
     if (ret  >= 0)
         *tid = v->rid, *beg = v->pos, *end = v->pos + v->rlen;
+    /*bcf data read may need conversion to vcf44 phasing format, as header is
+    not availble here, it has to be done after this returns the data*/
     return ret;
 }
 
@@ -2300,6 +2442,11 @@ 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
@@ -3225,7 +3372,7 @@ static int vcf_parse_format_fill5(kstring_t *s, const bcf_hdr_t *h, bcf1_t *v,
                         is_phased = (*t == '|');
                         if (*t != '|' && *t != '/') break;
                     }
-                    if (ver >= VCF44 && !phasingprfx) {
+                    if (!phasingprfx) { //get GT in v44 way when no prefixed phasing
                         /* no explicit phasing for 1st allele, set based on
                          other alleles and ploidy */
                         if (ploidy == 1) {  //implicitly phased