python pgenlib RuntimeError: read_alleles_range() error 6

[Kangcheng Hou]

Dear Chris,

We got the following message using the latest pgenlib read_alleles_range

 pgen.read_alleles_range(snp_start, snp_stop, geno)
  File "pgenlib.pyx", line 894, in pgenlib.PgenReader.read_alleles_range
  File "pgenlib.pyx", line 922, in pgenlib.PgenReader.read_alleles_range
RuntimeError: variant_idx 14875 read_alleles_range() error 6

The PLINK2 file was converted from VCF (obtained from michigan imputation server) with 

  --maf 0.005
  --make-pgen
  --max-alleles 2
  --out ${out}
  --rm-dup exclude-all
  --snps-only
  --vcf ${vcf}

We checked the .pvar file for 14875 ID and it looks normal.

My observations are:

1. PLINK2 binary is working properly with this file.

2. pgen.read_alleles_range was applied iteratively across multiple blocks in my application, and this function is working fine for other block of this file. So it is a problem to this variant.

3. I have applied the same code to other PLINK2 files and they were all working fine.

Would be great to hear your thoughts on this.

Best,

Kangcheng

[Christopher Chang]

Can you send me a dataset I can replicate this error with?

[Kangcheng Hou]

Hi Chris,

We are able to resolve this with different computing environment. This seems to be a platform related thing so may be harder to replicate.

This is just for your information and we are happy as is now.

Kangcheng

[Kangcheng Hou]

Hi Chris,

I am resending a replicable example just in case you didn't receive the one I last sent.

Steps to replicate: unzip the folder and run "python replicate.py". I am able to replicate read_alleles_range() error 6 in both OS X and linux.

Let me know anything i can help with this.

Kangcheng

[Christopher Chang]

Ok, the issue is that "hardcall_phase_present = True" was missing from the PgenWriter initialization statement.  I will look into making this error obvious.

[Christopher Chang]

I have posted an updated python package that catches this error.

[Kangcheng Hou]

Thanks for your quick response.

I ran the following code using the same data

import numpy as np
import pgenlib

input_path = "plink2.pgen"
out_path = "plink2.out.pgen"

with pgenlib.PgenReader(bytes(input_path, "utf8")) as pgen:
    n_snp = pgen.get_variant_ct()
    n_indiv = pgen.get_raw_sample_ct()
    print(f"Reading: {n_snp}, {n_indiv}")

    geno = np.empty([n_snp, pgen.get_raw_sample_ct() * 2], dtype=np.int32)
    pgen.read_alleles_range(0, n_snp, geno)

with pgenlib.PgenWriter(
    out_path.encode("utf-8"),
    sample_ct=n_indiv,
    variant_ct=n_snp,
    nonref_flags=False,
    hardcall_phase_present=True,
) as writer:
    writer.append_alleles_batch(geno, all_phased=True)

with pgenlib.PgenReader(bytes(out_path, "utf8")) as pgen:
    n_snp = pgen.get_variant_ct()
    n_indiv = pgen.get_raw_sample_ct()
    print(f"Reading again: {n_snp}, {n_indiv}")
    geno = np.empty([n_snp, pgen.get_raw_sample_ct() * 2], dtype=np.int32)
    pgen.read_alleles_range(0, n_snp, geno)

and get

RuntimeError: append_alleles_batch called with all_phased True, but PgenWriter was constructed with hardcall_phase_present False

Did I miss anything? Thanks in advance.

- Kangcheng

--
You received this message because you are subscribed to a topic in the Google Groups "plink2-users" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/plink2-users/5Xxa-DDcXRY/unsubscribe.
To unsubscribe from this group and all its topics, send an email to plink2-users...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/plink2-users/af4291dd-3302-4f66-8fad-8f7b0e4dd7bbn%40googlegroups.com.

[Christopher Chang]

No, I forgot to propagate the flag to the python object, sorry about that.  Posted another fix to GitHub.