OpenMendel / BGEN.jl

MIT License
10 stars 0 forks source link

Cannot get ref/alt allele label for haplotype data #15

Closed biona001 closed 3 years ago

biona001 commented 3 years ago

I think I'm running into the same problem as #7.

After parsing a variant's haplotypes, I want to additionally store ref/alt label (ie call major_allele and minor_allele), but I get Error: minor_allele_dosage!() must be called before minor_allele(). But calling minor_allele_dosage!() throws another error.

MWE

using BGEN
b = Bgen(BGEN.datadir("haplotypes.bgen"))
variants = parse_variants(b; from_bgen_start=true)
v = variants[1]

dose = probabilities!(b, v) # I will parse `dose` using some other routine to get haplotypes

julia> major_allele(v) # now try to get REF allele
┌ Error: `minor_allele_dosage!()` must be called before `minor_allele()`
└ @ BGEN ~/.julia/dev/BGEN/src/variant.jl:92

Try calling minor_allele_dosage!() throws not phased error

julia> minor_allele_dosage!(b, v)
ERROR: AssertionError: p.phased == 0
Stacktrace:
 [1] ref_allele_dosage!(b::Bgen, v::Variant; T::Type, mean_impute::Bool, clear_decompressed::Bool, data::Nothing, decompressed::Nothing, is_decompressed::Bool)
   @ BGEN ~/.julia/dev/BGEN/src/genotypes.jl:546
 [2] minor_allele_dosage!(b::Bgen, v::Variant; T::Type, mean_impute::Bool, clear_decompressed::Bool, data::Nothing, decompressed::Nothing, is_decompressed::Bool)
   @ BGEN ~/.julia/dev/BGEN/src/genotypes.jl:603
 [3] minor_allele_dosage!(b::Bgen, v::Variant)
   @ BGEN ~/.julia/dev/BGEN/src/genotypes.jl:589
 [4] top-level scope
   @ REPL[57]:1

but variant is phased:

julia> phased(v) == true
true

Here is another test data (VCF haplotypes converted to BGEN using qctools v2) Archive.zip

kose-y commented 3 years ago

Resolved by 8c2d462ea31e914d3b74319361a5eeac272df59f. Now minor_allele_dosage!() and ref_allele_dosage!() also work with biallelic phased data.

biona001 commented 3 years ago

thanks!!