Ragged haplotype sequences are unexpectedly truncated at the 3' end

[VincentLHH]

Description
I am encountering an issue where the haplotype sequences generated by genvarloader (v0.22.0) appear to be truncated at the end compared to the reference sequence length.

While the variants (indels) within the sequence are correctly applied, the resulting strings are missing bases at the 3' end. This results in a sequence length that does not match the theoretically calculated length based on the VCF indels.

Environment

• genvarloader version: 0.22.0
• Reference: hg38

Code to Reproduce

import genvarloader as gvl
from tqdm import tqdm
import json

DEBUG = True
OUTPUT_FILE = "output_haplotypes.json"

# 1. Create GVL file
gvl_path = "test_APOB.gvl"
gvl.write(
    path=gvl_path,
    bed="test.bed",
    variants="test.vcf.gz",
)

# 2. Open dataset and fetch ragged haplotypes
dataset = gvl.Dataset.open("test_APOB.gvl", reference="/ref/hg38.fa")
dataset = dataset.with_len("ragged").with_seqs("haplotypes")

genes = dataset.regions["name"].to_list()
samples = dataset.samples

haplotypes_dict = {}

for sample_idx, sample_name in enumerate(tqdm(samples, desc="Processing samples")):
    haplotypes_dict[sample_name] = {}
    for region_idx, gene_name in enumerate(genes):
        haps = dataset[region_idx, sample_idx]
        seq1 = haps[0].decode()
        seq2 = haps[1].decode()
        
        haplotypes_dict[sample_name][gene_name] = [seq1, seq2]

        if DEBUG:
            tqdm.write(f"{sample_name} _ {gene_name}: Strand1 len {len(seq1)}, Strand2 len {len(seq2)}")

with open(OUTPUT_FILE, "w") as f:
    json.dump(haplotypes_dict, f, indent=2)

Observed Behavior
The generated sequences are consistently shorter than the reference sequence at the end.

• Reference Sequence Length: 42645 bp
• Haplotype 1 Length: 42636 bp (9 bp shorter than ref)
• Haplotype 2 Length: 42645 bp (Matches ref length, but content is truncated)

Expected Behavior
Based on the variants in the phased VCF (specifically the indels), I calculated the expected lengths to be:

• Expected Haplotype 1: 42641 bp
• Expected Haplotype 2: 42647 bp

Furthermore, upon manual inspection of the sequence tails, bases are missing compared to the reference:

Reference Tail: ...GCAGGTCCCGGTGGGAAT
Haplotype 1 Tail: ...GCAGGTCCCGGTG (5 nt loss)
Haplotype 2 Tail: ...GCAGGTCCCGGTGGGA (2 nt loss)

The variants in the upstream region match perfectly, but the truncation occurs randomly at the very end of the string.

Question
Is this a bug regarding the "ragged" length setting, or is there a specific parameter I need to set to ensure the full sequence context is retrieved?

---

[d-laub]

Hi Vincent, thanks for opening this issue. GVL is continuously tested against bcftools consensus for correctness, however it's possible there are edge cases not covered. How did you compute the expected haplotype lengths? Can you confirm whether the variants passed to gvl.write were normalized and atomized? i.e. they were passed through bcftools norm -a --atom-overlaps . -m -any --multi-overlaps . -f /ref/hg38.fa Also, are you able to share test.bed and test.vcf.gz? Thanks!

[VincentLHH]

Thank you for the detailed guidance.

I have run the full bcftools norm command you suggested, but the issue persists. Below is a summary of my test setup and findings.

1. Test BED (0‑based)

$ cat test.bed
chr2    21001428        21044073        APOB

Region length = 42645 bp.

2. Test VCF

The VCF was extracted from a larger file using bcftools view -R test.bed -s sample_name.
For privacy I cannot share the full VCF, but I can show all variants that affect the two haplotypes for this sample.

All non‑reference variants (both haplotypes):

$ bcftools view -H test.vcf.gz | grep -v '0|0'

chr2    21008652        chr2_21008652_G_A       G       A       .       PASS    AC=2;AN=2       GT      1|1
chr2    21009931        chr2_21009931_T_C       T       C       .       PASS    AC=2;AN=2       GT:PP   1|1:.
chr2    21013610        chr2_21013610_A_G       A       G       .       PASS    AC=2;AN=2       GT:PP   1|1:.
chr2    21014672        chr2_21014672_G_A       G       A       .       PASS    AC=2;AN=2       GT      1|1
chr2    21016025        chr2_21016025_T_C       T       C       .       PASS    AC=2;AN=2       GT      1|1
chr2    21017492        chr2_21017492_A_G       A       G       .       PASS    AC=1;AN=2       GT      1|0
chr2    21018582        chr2_21018582_A_T       A       T       .       PASS    AC=2;AN=2       GT:PP   1|1:.
chr2    21018633        chr2_21018633_A_G       A       G       .       PASS    AC=1;AN=2       GT      1|0
chr2    21022457        chr2_21022457_G_T       G       T       .       PASS    AC=1;AN=2       GT      1|0
chr2    21022495        chr2_21022495_A_G       A       G       .       PASS    AC=1;AN=2       GT      1|0
chr2    21023434        chr2_21023434_G_C       G       C       .       PASS    AC=1;AN=2       GT      1|0
chr2    21028042        chr2_21028042_G_A       G       A       .       PASS    AC=1;AN=2       GT      0|1
chr2    21028168        chr2_21028168_T_G       T       G       .       PASS    AC=1;AN=2       GT      0|1
chr2    21029187        chr2_21029187_C_T       C       T       .       PASS    AC=1;AN=2       GT:PP   0|1:1
chr2    21029493        chr2_21029493_GAGGA_G   GAGGA   G       .       PASS    AC=1;AN=2       GT      0|1
chr2    21029493        chr2_21029493_G_GAGGA   G       GAGGA   .       PASS    AC=1;AN=2       GT      1|0
chr2    21029496        chr2_21029496_GAA_*     GAA     *       .       PASS    AC=1;AN=2       GT      0|1
chr2    21032701        chr2_21032701_C_CCATCT  C       CCATCT  .       PASS    AC=1;AN=2       GT      0|1
chr2    21036740        chr2_21036740_A_G       A       G       .       PASS    AC=1;AN=2       GT      0|1
chr2    21037212        chr2_21037212_G_A       G       A       .       PASS    AC=1;AN=2       GT      0|1
chr2    21037656        chr2_21037656_G_GT      G       GT      .       PASS    AC=2;AN=2       GT:PP   1|1:.
chr2    21037729        chr2_21037729_G_A       G       A       .       PASS    AC=1;AN=2       GT      0|1
chr2    21038635        chr2_21038635_T_C       T       C       .       PASS    AC=1;AN=2       GT      0|1
chr2    21039254        chr2_21039254_A_T       A       T       .       PASS    AC=1;AN=2       GT      0|1
chr2    21039555        chr2_21039555_G_A       G       A       .       PASS    AC=1;AN=2       GT      0|1
chr2    21040767        chr2_21040767_T_G       T       G       .       PASS    AC=2;AN=2       GT      1|1
chr2    21041862        chr2_21041862_T_C       T       C       .       PASS    AC=1;AN=2       GT      1|0
chr2    21042065        chr2_21042065_T_G       T       G       .       PASS    AC=1;AN=2       GT      1|0
chr2    21042269        chr2_21042269_C_A       C       A       .       PASS    AC=2;AN=2       GT      1|1
chr2    21043292        chr2_21043292_C_T       C       T       .       PASS    AC=1;AN=2       GT      1|0
chr2    21043902        chr2_21043902_GGCAGCGCCA_G      GGCAGCGCCA      G       .       PASS    AC=1;AN=2       GT      1|0
chr2    21043906        chr2_21043906_GCGCCA_*  GCGCCA  *       .       PASS    AC=1;AN=2       GT      1|0

Variants on haplotype 1 (first allele non‑reference):

$ bcftools view -H test.vcf.gz | grep -v '0|0' | grep -v '0|1'

chr2    21008652        chr2_21008652_G_A       G       A       ... GT      1|1
chr2    21009931        chr2_21009931_T_C       T       C       ... GT:PP   1|1:.
chr2    21013610        chr2_21013610_A_G       A       G       ... GT:PP   1|1:.
chr2    21014672        chr2_21014672_G_A       G       A       ... GT      1|1
chr2    21016025        chr2_21016025_T_C       T       C       ... GT      1|1
chr2    21017492        chr2_21017492_A_G       A       G       ... GT      1|0
chr2    21018582        chr2_21018582_A_T       A       T       ... GT:PP   1|1:.
chr2    21018633        chr2_21018633_A_G       A       G       ... GT      1|0
chr2    21022457        chr2_21022457_G_T       G       T       ... GT      1|0
chr2    21022495        chr2_21022495_A_G       A       G       ... GT      1|0
chr2    21023434        chr2_21023434_G_C       G       C       ... GT      1|0
chr2    21029493        chr2_21029493_G_GAGGA   G       GAGGA   ... GT      1|0
chr2    21037656        chr2_21037656_G_GT      G       GT      ... GT:PP   1|1:.
chr2    21040767        chr2_21040767_T_G       T       G       ... GT      1|1
chr2    21041862        chr2_21041862_T_C       T       C       ... GT      1|0
chr2    21042065        chr2_21042065_T_G       T       G       ... GT      1|0
chr2    21042269        chr2_21042269_C_A       C       A       ... GT      1|1
chr2    21043292        chr2_21043292_C_T       C       T       ... GT      1|0
chr2    21043902        chr2_21043902_GGCAGCGCCA_G      GGCAGCGCCA      G       ... GT      1|0
chr2    21043906        chr2_21043906_GCGCCA_*  GCGCCA  *       ... GT      1|0

The indels that change length on haplotype 1 are:

• chr2 21029493 G→GAGGA (+4 bp)
• chr2 21037656 G→GT (+1 bp)
• chr2 21043902 GGCAGCGCCA→G (-9 bp)
• chr2 21043906 GCGCCA→* (-6bp but overlaps with the -9bp deletion above)

Thus the expected length for haplotype 1 = 42645 + 4 + 1 - 9 = 42641.

No variant is located near the region end (21044073) – the last variant is at 21043906 (167 bp upstream). So the 3′ end should be intact.

3. What I observed with GenVarLoader (v0.22.0)

Using ragged output (with_len("ragged")), the obtained lengths were:

• Haplotype 1: 42636 (expected 42641 → 5 bp too short)
• Haplotype 2: 42645 (expected 42647 → 2 bp too short)

The terminal sequences compared to the reference:

• Reference end: ...GCAGGTCCCGGTGGGAAT
• Haplotype 1 end: ...GCAGGTCCCGGTG (missing GGAAT)
• Haplotype 2 end: ...GCAGGTCCCGGTGGGA (missing AT)

All internal variants matched correctly; only the very 3′ end was truncated.

4. Preprocessing check

I ran the full bcftools norm command as you suggested:

bcftools norm -a --atom-overlaps . -m -any --multi-overlaps . -f hg38.fa test.vcf.gz -Oz -o test_norm.vcf.gz

The output:

Lines total/split/joined/realigned/mismatch_removed/dup_removed/skipped: 3894/0/0/0/0/0/0

This shows the VCF was already normalized and atomized.

I then rebuilt the GVL dataset using test_norm.vcf.gz and re‑extracted the haplotypes. The results were identical to those from the original VCF (diff haplotypes_test.json haplotypes_test_norm.json produced no output). The terminal truncation persists.

5. Conclusion

The VCF is properly phased, normalized, and atomized. The length calculations based on the VCF match the reference minus the exact indels. However, GenVarLoader produces haplotypes that are shorter than expected at the 3′ end, and the missing bases do not correspond to any recorded variant.

Could this be a bug in the ragged haplotype reconstruction when the last variant(s) are indels? Or is there an internal boundary‑trimming mechanism that I have missed? Any debugging suggestions would be greatly appreciated.

Thank you again for your time.