Phased GT encoding

[yuw444]

I am working a project using <htslib/vcf.h>. In particular, I am trying to access gene type (GT) info. I've shared my findings below.

1. GT is an encoded array with allele and phase info.
2. Idx is the index for alleles w.r.t the ref allele(0). For instance,

idx = 0(ref), 1, 2, 3(alt)

3. when phased(|), the corresponding encoded values are

val = 3, 5, 7, 9;

4. when unphased (), the corresponding encoded values are

val = 2, 4, 6, 8;

    // Macros for setting genotypes correctly, for use with bcf_update_genotypes only; idx corresponds
    // to VCF's GT (1-based index to ALT or 0 for the reference allele) and val is the opposite, obtained
    // from bcf_get_genotypes() below.
    #define bcf_gt_phased(idx)      (((idx)+1)<<1|1)
    #define bcf_gt_unphased(idx)    (((idx)+1)<<1)
    #define bcf_gt_missing          0
    #define bcf_gt_is_missing(val)  ((val)>>1 ? 0 : 1)
    #define bcf_gt_is_phased(idx)   ((idx)&1)   // ??? should it be  bcf_gt_is_phased(val)???
    #define bcf_gt_allele(val)      (((val)>>1)-1)

https://github.com/samtools/htslib/blob/c705bec2af9caedacc6ba80dfba7fe625eb5ba5c/htslib/vcf.h#L1019C1-L1028C1

Later, one can also decode its val back to idx and phase. However, when I am parsing an example phased bcf file as attached, I found the following

#CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  ID1       ID2
chr1    54490   .       G       A       .       .       GC=0.394015     GT:AD:Ldev:Bdev:AS      1|0:19,19:0:0:0 ./.:.:.:.:.
chr1    634553  .       G       A       .       .       GC=0.511222     GT:AD:Ldev:Bdev:AS      ./.:.:.:.:.     0|1:14,21:0.002:0.158:-1

1. The first record chr1:54490, GT=1|0; ./., according to the above encode algorithm, my projected values should be
   (((1) + 1) << 1 | 1) = 2 * 2 + 1 = 5,
(((0) + 1) << 1 | 1) = 1 * 2 + 1 = 3.
   However, the observed values are 4 and 3 from (int8_t *)fmt->p.

2. The second record chr1:634553, GT=./.;0:1. Similarly, the projected values are 3 and 5. However, the observed values are 2 and 5.

This leads to my conjecture that (int8_t *)fmt->p [0] is not been considered when encoding somewhere.

Then, I checked for the unphased same bcf file, the encoded values are (4, 2) and (2,4), which is consistent with the above code snippets.

#CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  ID1       ID2
chr1    54490   .       G       A       .       .       GC=0.394015     GT:AD:Ldev:Bdev:AS      1/0:19,19:0:0:0 ./.:.:.:.:.
chr1    634553  .       G       A       .       .       GC=0.511222     GT:AD:Ldev:Bdev:AS      ./.:.:.:.:.     0/1:14,21:0.002:0.158:-1

Albeit, the current coding scheme doesn’t affect bcf_gt_allele at all, and has some effect when evaluating bcf_gt_is_phased(val).

I’m not sure if there is any special consideration for position 0. Otherwise, I would be more than happy to work with you to look into the details. I have attached the code I used to read.

test_read_bcf.txt

[daviesrob]

Your code looks incorrect. See bcf_format_gt() to see how to read GT values (up to VCF4.3). Compared with your code, note that only one value is accessed for each GT. It is also not valid C to cast the pointer into the GT data to a wider type, as it may not be correctly aligned. That's why bcf_format_gt() internally uses uint8_t pointers and functions like le_to_i16 or le_to_i32 to access the data.

[yuw444]

@daviesrob , please see this bcf_all_phased. I feel I used the same approach as what it does. However, for loop can still assign sample_phased = 1 as long as p[1] is odd even if p[0] is not odd.

[daviesrob]

Ah, I see your code was a bit specialised to having two genotypes per sample. If you look at the code for bcf_all_phased that you linked, you'll see it goes through all the samples:

        for (isample=0; isample<line->n_sample; isample++)

and then all the genotypes for that sample in an inner loop:

                for (i=0; i<fmt_ptr->n; i++) \

as VCF supports samples with ploidy greater than two.

Apart from assuming a diploid record, your code is actually working as expected. The first genotype always appearing unphased comes from the VCF encoding, where the phasing marker always appears between the allele numbers: 0/1/2, 0|1/2. As the number of phasing markers is always one less than the number of genotypes for the sample, the first is currently always stored as unphased. Note that this will change in a future version of HTSlib, for VCFv4.4 and later files. VCFv4.4 extends the definition of GT to allow an optional leading phasing indicator, so records like |0/1/2 will become possible. So once #1861 has been merged, you will see the values you originally expected, although only for VCFv4.4 or later files.

[yuw444]

I see your point now. As in the current version, the first position will always be unphased despite the phased info for the rest.

Looking forward to the new version.