TypeError when running SeqrMTToESTask

[tomiles]

When running python3 -m seqr_loading SeqrMTToESTask --local-scheduler I run into an error.
Looking at the source code I can't figure out why this occurs. Hopefully someone here can provide some assistance on how I can further debug this.

Version info:

• Hail version 0.2.94-f0b38d6c436f
• Python 3.7.1

Error log:

ERROR: [pid 6174] Worker Worker(salt=664600238, workers=1, host=3da1e4db6e23, username=root, pid=6174) failed    SeqrVCFToMTTask(source_paths=/input_vcfs/WGS/Callset-gatk-haplotype-joint-annotated-filtered.vcf.gz, dest_path=/input_vcfs/WGS/Callset-gatk-haplotype-joint-annotated-filtered.mt, genome_version=38, vep_runner=VEP, reference_ht_path=/seqr-reference-data/GRCh38/combined_reference_data_grch38.ht, clinvar_ht_path=/seqr-reference-data/GRCh38/clinvar.GRCh38.ht, hgmd_ht_path=None, sample_type=WGS, dont_validate=False, dataset_type=VARIANTS, remap_path=, subset_path=, vep_config_json_path=/vep_configs/vep-GRCh38-loftee.json, grch38_to_grch37_ref_chain=gs://hail-common/references/grch38_to_grch37.over.chain.gz)
Traceback (most recent call last):
  File "/usr/local/lib/python3.7/site-packages/luigi/worker.py", line 199, in run
    new_deps = self._run_get_new_deps()
  File "/usr/local/lib/python3.7/site-packages/luigi/worker.py", line 141, in _run_get_new_deps
    task_gen = self.task.run()
  File "/hail-elasticsearch-pipelines/luigi_pipeline/seqr_loading.py", line 86, in run
    self.read_vcf_write_mt()
  File "/hail-elasticsearch-pipelines/luigi_pipeline/seqr_loading.py", line 110, in read_vcf_write_mt
    overwrite=True).select_annotated_mt()
  File "/hail-elasticsearch-pipelines/luigi_pipeline/lib/model/base_mt_schema.py", line 162, in annotate_all
    func_ret = annotation.fn(self)
  File "/hail-elasticsearch-pipelines/luigi_pipeline/lib/model/seqr_mt_schema.py", line 220, in ac
    return self.mt.info.AC[self.mt.a_index-1]
TypeError: 'Int32Expression' object is not subscriptable

[hanars]

The command you include here is not the recommended way to run the seqr pipeline. Please make sure you are running the latest docker image, have followed all the set up instructions for pipeline reference data, and are running one of the supported loading commands as outlined in the latest version of our instructions: https://github.com/broadinstitute/seqr/blob/master/deploy/LOCAL_INSTALL.md#annotating-and-loading-vcf-callsets

[tomiles]

I am running according to the docs you reference, using the load_data.sh script in the latest docker container but that just wraps python3 -m seqr_loading SeqrMTToESTask --local-scheduler ... So that is why that command is in my post.

[hanars]

I see. In the future, including the command you actually run is much more helpful than putting in part of the command that is executed. Please provide the actual command you ran, including all the arguments, and the full logs.

Without the full logs I can't do much more, but this might be an issue with the input VCF format. Can you confirm that you are running the pipeline on a GATK joint called VCF?

[tomiles]

The vcf is indeed a GATK joint called VCF. I could not find anything relevant in the full log besides the selected part I posted above. Is there a way to enable more detailed logging, so we can pinpoint for which variant in the vcf the AC field triggers this Type error in seqr_mt_schema.py? (besides editing the python source to add some debug print statements)

[tomiles]

FYI the TypeError: 'Int32Expression' object is not subscriptable was triggered because we tried to load a vcf that was decomposed (usingvt) and that has AC/AN/AF field in the header lines as Number=1 while hail seems to expectNumber=A for those fields.

[hanars]

Please provide the actual command you ran, including all the arguments, and the full logs. Logs are available by running docker-compose logs pipeline-runner. Without this information I can not help you

[tomiles]

Since there was nothing informative in the logs I reran the analyses with a newly "clean" called vcf (without existing annotations). But now we hit another error:

ERROR: [pid 70616] Worker Worker(salt=594408665, workers=1, host=608cf974c570, username=root, pid=70616) failed    SeqrVCFToMTTask(source_paths=/input_vcfs/WGS/Anonymous-joint-filtered_debugline.vcf.gz, dest_path=/input_vcfs/WGS/Anonymous-joint-filtered_debugline.mt, genome_version=38, vep_runner=VEP, reference_ht_path=/seqr-reference-data/GRCh38/combined_reference_data_grch38.ht, clinvar_ht_path=/seqr-reference-data/GRCh38/clinvar.GRCh38.ht, hgmd_ht_path=None, sample_type=WGS, dont_validate=True, dataset_type=VARIANTS, remap_path=, subset_path=, vep_config_json_path=/vep_configs/vep-GRCh38-loftee.json, grch38_to_grch37_ref_chain=gs://hail-common/references/grch38_to_grch37.over.chain.gz)
Traceback (most recent call last):
  File "/usr/local/lib/python3.7/site-packages/luigi/worker.py", line 199, in run
    new_deps = self._run_get_new_deps()
  File "/usr/local/lib/python3.7/site-packages/luigi/worker.py", line 141, in _run_get_new_deps
    task_gen = self.task.run()
  File "/hail-elasticsearch-pipelines/luigi_pipeline/seqr_loading.py", line 86, in run
    self.read_vcf_write_mt()
  File "/hail-elasticsearch-pipelines/luigi_pipeline/seqr_loading.py", line 118, in read_vcf_write_mt
    mt.write(self.output().path, stage_locally=True, overwrite=True)
  File "<decorator-gen-1154>", line 2, in write
  File "/usr/local/lib/python3.7/site-packages/hail/typecheck/check.py", line 577, in wrapper
    return __original_func(*args_, **kwargs_)
  File "/usr/local/lib/python3.7/site-packages/hail/matrixtable.py", line 2556, in write
    Env.backend().execute(ir.MatrixWrite(self._mir, writer))
  File "/usr/local/lib/python3.7/site-packages/hail/backend/py4j_backend.py", line 104, in execute
    self._handle_fatal_error_from_backend(e, ir)
  File "/usr/local/lib/python3.7/site-packages/hail/backend/backend.py", line 189, in _handle_fatal_error_from_backend
    raise HailUserError(message_and_trace) from None
hail.utils.java.HailUserError: Error summary: HailException: array index out of bounds: index=2, length=2

When I narrow this down to where in the vcf this exception gets triggered it seems to only happen for variants on chrX, chrY and chrM. Loading the vcf without those 3 contigs works perfectly (chr1 -> 22).

Test VCF data:

#CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  SAMPLE_M       SAMPLE_F

Working vcf line:

chr1    10237   rs1490526935    A       C       20      PASS    AC=2;AF=1;AN=2;DB;DP=46;ExcessHet=3.0103;FS=0;MLEAC=2;MLEAF=1;MQ=9.81;MQ0=0;QD=20;SOR=1.609     GT:AD:DP:GQ:PGT:PID:PL  1/1:0,1:1:3:1|1:10237_A_C:45,3,0        ./.:40,0:40:.:.:.:0,0,0

Failing vcf lines:

chrY    2791306 rs111319867     A       C       364     PASS    AC=1;AF=1;AN=1;DB;DP=11;FS=0;MLEAC=1;MLEAF=1;MQ=60;MQ0=0;QD=33.09;SOR=0.859     GT:AD:DP:GQ:PL  1:0,11:11:99:394,0      .:.:.:.:.
chrM    73      rs869183622     A       G       63717.8 PASS    AC=1;AF=0.5;AN=2;DB;DP=2122;FS=0;MLEAC=1;MLEAF=0.5;MQ=57.42;MQ0=0;QD=34.24;SOR=0.922    GT:AD:DP:GQ:PL  0:502,0:502:99:0,1800   1:0,1604:1604:99:63749,0

Any idea why this HailException would get triggered only on chrX, chrY and chrM variants in the vcf? could this be a ploidy handling issue?

[hanars]

@williamphu have you ever encountered this issue before, or know why it might be happening?

[tomiles]

included here a log file when test loading with an artificially created vcf only including 1 chrX variant triggering the array index out of bounds :

hail-20220602-1650-0.2.94-f0b38d6c436f.log

[tomiles]

Also tried with different real vcf's that trigger this error and filtered out chrX, chrY and chrM variants beforehand with bcftools and this makes them import perfectly fine.

[williamphu]

I've run into this error before, and it is related to ploidy - for context, the issue is that our loading pipeline uses Hail's function for splitting multi-allelic sites: https://hail.is/docs/0.2/methods/genetics.html#hail.methods.split_multi_hts - which expects three scores for PL (and might expect a specific number for other scores, I've run into the PL issue first). For variants with a haploid call, our pipeline will throw an error because it's expecting scores that don't exist for a haploid call. Currently we don't have a good fix for this other than not including haploid calls - digging into the source for the function https://hail.is/docs/0.2/_modules/hail/methods/statgen.html#split_multi_hts, you might be able to work around this by removing the PL field from the VCF - but our downstream code might be expecting a PL score.

[williamphu]

To confirm this, can you try filtering out your haploid calls in your X, Y, and M variants?

[tomiles]

You mean to retry with the PL field removed? I will edit the file and rerun.

[tomiles]

To confirm removing the PL field does indeed make the pipeline process the variants successfully.

Is there some known best practice for GATK Haplotype Caller to produce a more compatible vcf for these haploid calls and avoid having to do the extra bcftool tweaking/filtering ? Will have to discuss this with the scientist doing the variant interpretation if they can tolerate loosing the PL info for variants on those chromosomes.

[hanars]

For what its worth the pipeline later drops the PL field entirely so it was never going to be shown in seqr anyways. So even if you left it in the seqr users would not have see it

[tomiles]

Good to know, then I'll just drop the field on the whole vcf before importing.