"found allele outside of expected range" error

minku · November 30, 2020, 2:27am

Hi there!

whenever I try to mt.write(), this error shows up:

HailException: found allele outside of expected range [0, 2]: 3

My matrix table has undergone: split_multi_hts(), vep(), and variant_qc().
I had no problem show()ing this matrixtable, but only when I try to write() onto my disk, this error appears.

I want to know:

What does the error message mean? Is it saying that although the mt.alleles were composed of 0, 1, and 2 alleles → the mt.GT is showing genotype that contains allele that is out of this range?(such as 3?)
What should I do to fix this?

tpoterba · November 30, 2020, 2:57pm

Lazy evaluation

One of the most challenging pieces of Hail’s learning curve is its lazy evaluation. This means that Hail doesn’t actually run any computation until it hits a certain kind of method – write, aggregate, export, show.

In the following pipeline:

mt = hl.import_vcf(path)
mt.write('...')

If the VCF at path is malformed outside the header, you’d see an error in the second line, not the first, because Hail hasn’t actually parsed the VCF until the second line. This is not great from a usability perspective, but it affords us two advantages:

Scalability: Hail can process more data than can fit in memory
Performance: Hail can employ a query optimizer that rewrites full queries to do minimal work.

This particular error

I think this error in particular indicates that there is an unexpected allele somewhere, probably in split_multi_hts. Despite the lazy evaluation, we try to capture the Python stack trace of the original call site that triggered the error – can you paste the full stack trace (especially the python bit if there is one)?

minku · December 3, 2020, 1:38pm

Hi Tim!
Thank you so much for the clear explanation! I think I understand Hail better now.

I copy & pasted the stack trace below. Should I also put down the java stack trace?

---------------------------------------------------------------------------

   Traceback (most recent call last)
   /usr/local/lib/python3.6/dist-packages/IPython/core/formatters.py in __call__(self, obj)
       700                 type_pprinters=self.type_printers,
       701                 deferred_pprinters=self.deferred_printers)
   --> 702             printer.pretty(obj)
       703             printer.flush()
       704             return stream.getvalue()

   /usr/local/lib/python3.6/dist-packages/IPython/lib/pretty.py in pretty(self, obj)
       392                         if cls is not object \
       393                                 and callable(cls.__dict__.get('__repr__')):
   --> 394                             return _repr_pprint(obj, self, cycle)
       395 
       396             return _default_pprint(obj, self, cycle)

   /usr/local/lib/python3.6/dist-packages/IPython/lib/pretty.py in _repr_pprint(obj, p, cycle)
       698     """A pprint that just redirects to the normal repr function."""
       699     # Find newlines and replace them with p.break_()
   --> 700     output = repr(obj)
       701     lines = output.splitlines()
       702     with p.group():

   ~/.local/lib/python3.6/site-packages/hail/table.py in __repr__(self)
      1295 
      1296         def __repr__(self):
   -> 1297             return self.__str__()
      1298 
      1299         def data(self):

   ~/.local/lib/python3.6/site-packages/hail/table.py in __str__(self)
      1292 
      1293         def __str__(self):
   -> 1294             return self._ascii_str()
      1295 
      1296         def __repr__(self):

   ~/.local/lib/python3.6/site-packages/hail/table.py in _ascii_str(self)
      1318                 return s
      1319 
   -> 1320             rows, has_more, dtype = self.data()
      1321             fields = list(dtype)
      1322             trunc_fields = [trunc(f) for f in fields]

   ~/.local/lib/python3.6/site-packages/hail/table.py in data(self)
      1302                 row_dtype = t.row.dtype
      1303                 t = t.select(**{k: hl._showstr(v) for (k, v) in t.row.items()})
   -> 1304                 rows, has_more = t._take_n(self.n)
      1305                 self._data = (rows, has_more, row_dtype)
      1306             return self._data

   ~/.local/lib/python3.6/site-packages/hail/table.py in _take_n(self, n)
      1449             has_more = False
      1450         else:
   -> 1451             rows = self.take(n + 1)
      1452             has_more = len(rows) > n
      1453             rows = rows[:n]

   <decorator-gen-1119> in take(self, n, _localize)

   ~/.local/lib/python3.6/site-packages/hail/typecheck/check.py in wrapper(__original_func, *args, **kwargs)
       612     def wrapper(__original_func, *args, **kwargs):
       613         args_, kwargs_ = check_all(__original_func, args, kwargs, checkers, is_method=is_method)
   --> 614         return __original_func(*args_, **kwargs_)
       615 
       616     return wrapper

   ~/.local/lib/python3.6/site-packages/hail/table.py in take(self, n, _localize)
      2119         #""" 
      2120 
   -> 2121         return self.head(n).collect(_localize)
      2122 
      2123     @typecheck_method(n=int)

   <decorator-gen-1113> in collect(self, _localize)

   ~/.local/lib/python3.6/site-packages/hail/typecheck/check.py in wrapper(__original_func, *args, **kwargs)
       612     def wrapper(__original_func, *args, **kwargs):
       613         args_, kwargs_ = check_all(__original_func, args, kwargs, checkers, is_method=is_method)
   --> 614         return __original_func(*args_, **kwargs_)
       615 
       616     return wrapper

   ~/.local/lib/python3.6/site-packages/hail/table.py in collect(self, _localize)
      1918         e = construct_expr(rows_ir, hl.tarray(t.row.dtype))
      1919         if _localize:
   -> 1920             return Env.backend().execute(e._ir)
      1921         else:
      1922             return e

   ~/.local/lib/python3.6/site-packages/hail/backend/py4j_backend.py in execute(self, ir, timed)
        96                 raise HailUserError(message_and_trace) from None
        97 
   ---> 98             raise e

   ~/.local/lib/python3.6/site-packages/hail/backend/py4j_backend.py in execute(self, ir, timed)
        72         # print(self._hail_package.expr.ir.Pretty.apply(jir, True, -1))
        73         try:
   ---> 74             result = json.loads(self._jhc.backend().executeJSON(jir))
        75             value = ir.typ._from_json(result['value'])
        76             timings = result['timings']

   ~/.local/lib/python3.6/site-packages/py4j/java_gateway.py in __call__(self, *args)
      1255         answer = self.gateway_client.send_command(command)
      1256         return_value = get_return_value(
   -> 1257             answer, self.gateway_client, self.target_id, self.name)
      1258 
      1259         for temp_arg in temp_args:

   ~/.local/lib/python3.6/site-packages/hail/backend/py4j_backend.py in deco(*args, **kwargs)
        30                 raise FatalError('%s\n\nJava stack trace:\n%s\n'
        31                                  'Hail version: %s\n'
   ---> 32                                  'Error summary: %s' % (deepest, full, hail.__version__, deepest), error_id) from None
        33         except pyspark.sql.utils.CapturedException as e:
        34             raise FatalError('%s\n\nJava stack trace:\n%s\n'

   FatalError: HailException: found allele outside of expected range [0, 2]: 3

hdjc90 · August 3, 2022, 9:01am

Hi @minku and @tpoterba did you ever find a solution to this problem? Many thanks

tpoterba · August 3, 2022, 11:47am

The issue here is that the alleles passed to the call_stats aggregator (typically mt.alleles, especially if this is coming from inside variant_qc) don’t match the alleles in a genotype call. This can happen if you have one ref/alt (2 alleles total) but have a 0/2 call (the 2 allele is out of range of the defined alleles).

hdjc90 · August 3, 2022, 1:09pm

Thanks @tpoterba i think that makes sense. What can cause this to happen? And Is there a workaround for this? Thanks in advance for your help

hdjc90 · August 3, 2022, 1:24pm

If it’s any help the variant caller used was illumina’s DRAGEN platform

tpoterba · August 3, 2022, 1:29pm

Here’s a way for you do figure out where the problem is –

mt = mt.annotate_rows(n_alleles = hl.len(mt.alleles))
print(mt.aggregate_entries(
  hl.agg.filter(hl.range(mt.GT.ploidy).any(lambda i: mt.GT[i] >= mt.n_alleles),
    hl.agg.take(hl.struct(locus=mt.locus, alleles=mt.alleles, sample=mt.s, entry=mt.entry)), 10)))

hdjc90 · August 3, 2022, 1:32pm

Thank you @tpoterba I will do this now. Much appreciated

hdjc90 · August 3, 2022, 1:54pm

tpoterba:

mt = mt.annotate_rows(n_alleles = hl.len(mt.alleles))
print(mt.aggregate_entries(
  hl.agg.filter(hl.range(mt.GT.ploidy).any(lambda i: mt.GT[i] >= mt.n_alleles),
    hl.agg.take(hl.struct(locus=mt.locus, alleles=mt.alleles, sample=mt.s, entry=mt.entry)), 10)))

Sorry if I have missed something here, but when I run this I get the following attribute error:

AttributeError: Table instance has no field, method, or property 'annotate_rows'
    Did you mean:
        Table methods: 'annotate', 'annotate_globals'

tpoterba · August 3, 2022, 1:55pm

Can you share the failing script? What are you calling call_stats or variant_qc on?

hdjc90 · August 3, 2022, 2:02pm

Of course:

import hail as hl

hl.init(tmp_dir='/...', spark_conf={

    "spark.executor.extraClassPath": "/usr/lib64/libblas.so:/usr/lib64/liblapack.so:/usr/lib64/libgslcblas.so:/usr/lib64/liblapack.so:/usr/lib64/libgslcblas.so:/usr/lib64/atlas/libsatlas.so.3"})

from hail import *

from pprint import pprint

import hail.expr.functions

import bokeh

from bokeh.io import show, output_notebook

from bokeh.layouts import gridplot

from bokeh.io import export_png

from bokeh.plotting import figure, output_file, save

mt = hl.import_vcf('/.../WGS_multisample.vcf.gz', reference_genome='GRCh38', force_bgz=True, skip_invalid_loci=True, array_elements_required=False)

mt = hl.split_multi(mt)
mt = hl.variant_qc(mt)
pedigree = hl.Pedigree.read('/.../trios.fam')
mt = hl.de_novo(mt, pedigree, mt.variant_qc.AF[1])

mt = mt.annotate_rows(n_alleles = hl.len(mt.alleles))
print(mt.aggregate_entries(
  hl.agg.filter(hl.range(mt.GT.ploidy).any(lambda i: mt.GT[i] >= mt.n_alleles),
    hl.agg.take(hl.struct(locus=mt.locus, alleles=mt.alleles, sample=mt.s, entry=mt.entry)), 10)))

Thanks again

tpoterba · August 3, 2022, 2:24pm

And if you put my two lines at the bottom there you see the “Table instance has no field, method, or property 'annotate_rows” error?

tpoterba · August 3, 2022, 2:24pm

oops, sorry, I didn’t scroll down…

tpoterba · August 3, 2022, 2:25pm

Ugh, the problem here is one of documentation, really. Instead of using hl.split_multi, you should use hl.split_multi_hts. The difference is that split_multi_hts updates genotype fields (GT/AD/PL/etc) but split_multi does not and leaves you to do this manually.

hdjc90 · August 3, 2022, 2:32pm

I see, thank you.

I have used hl.split_multi_hts previously but run into a different error that I wasn’t able to sort out either. Which led me to used hl.split_multi

This error is:

HailUserError: Error summary: HailException: array index out of bounds: index=3, length=3
------------
Hail stack trace:
  File "<ipython-input-3-fc98888815ec>", line 1, in <module>
    mt = hl.split_multi_hts(mt)

  File "/apps/genomics/hail/0.2/el7/AVX512/gnu-7.3/hail/methods/statgen.py", line 2322, in split_multi_hts
    (hl.range(0, 3).map(lambda i:

  File "/apps/genomics/hail/0.2/el7/AVX512/gnu-7.3/hail/methods/statgen.py", line 2326, in <lambda>
    ).map(lambda j: split.PL[j]))))))

  File "/apps/genomics/hail/0.2/el7/AVX512/gnu-7.3/hail/methods/statgen.py", line 2326, in <lambda>
    ).map(lambda j: split.PL[j]))))))

  File "/apps/genomics/hail/0.2/el7/AVX512/gnu-7.3/hail/expr/expressions/typed_expressions.py", line 481, in __getitem__
    return self._method("indexArray", self.dtype.element_type, item)

  File "/apps/genomics/hail/0.2/el7/AVX512/gnu-7.3/hail/expr/expressions/base_expression.py", line 596, in _method
    x = ir.Apply(name, ret_type, self._ir, *(a._ir for a in args))

  File "/apps/genomics/hail/0.2/el7/AVX512/gnu-7.3/hail/ir/ir.py", line 2138, in __init__
    self.save_error_info()

Originally I thought this was the inclusion of the <NON_REF> field shown below:

However when this is removed i still run into this issue.

Do you have any means to try and idenitfy this source of the error?

tpoterba · August 4, 2022, 11:49am

Could you show us what some examples of PLs look like for chr1:10002? Do they have 3 or 6 values?

hdjc90 · August 4, 2022, 12:14pm

They have 3 values for Chr1:10002 they appear like this :0,0,0:

hdjc90 · August 4, 2022, 12:26pm

for others like Chr1:10007 they are :0,8,816:

tpoterba · August 4, 2022, 12:44pm

Okay, looks like the G-indexing does not include the NON-REF allele, as expected.

Do you mean that you have a VCF without the NON_REF allele and it still has the same problem?

If the VCF still has the NON_REF allele included, let’s try removing it on the fly as the VCF is imported using the filter_replace arg:

mt = hl.import_vcf('/.../WGS_multisample.vcf.gz', reference_genome='GRCh38', force_bgz=True, skip_invalid_loci=True, array_elements_required=False, find_replace=(',<NON_REF>', ''))

mt = hl.split_multi_hts(mt)
mt = hl.variant_qc(mt)
mt.rows()._force_count() # forces evaluation to make sure the above runs without error

Topic		Replies	Views
Connection Errors Hail Query & hailctl	8	396	August 2, 2023
PL of haploid call induced error when split_multi_hts(), in VCF produced by Dragen pipeline Hail Query & hailctl	2	345	February 10, 2023
Write Bgen Fatal Error Hail Query & hailctl	3	336	July 17, 2023
Error index out of bounds Hail Query & hailctl	6	1085	July 27, 2022
VCF4.1, is.hail.utils.HailException: array index out of bounds Hail Query & hailctl	2	796	January 2, 2020

"found allele outside of expected range" error

Lazy evaluation

This particular error

Related Topics