genmod models can give a `KeyError` with `--phased`parameter if some variants lack region annotation

[sofiademmou]

genmod models gives a KeyError when used with multiple samples and --phased flag

Command

 genmod \
    models \
    --whole_gene --temp_dir ./ --phased  \
    --reduced_penetrance reduced_penetrance.tsv \
    --family_file FAM.ped \
    --processes 2 \
    --outfile test_data_chunk.00002.bed_snv_genmod_models_models.vcf \
    test_data_chunk.00002.bed_snv_genmod_annotate_annotate.vcf

Files used

FAM.ped.txt
reduced_penetrance.tsv
test_data_chunk.00002.bed_snv_genmod_annotate_annotate.vcf.txt

Error

> Traceback (most recent call last): 
 > File "/opt/conda/lib/python3.11/multiprocessing/process.py", line 314, in _bootstrap 
 > self.run() 
 > File "/opt/conda/lib/python3.11/site-packages/genmod/annotate_models/variant_annotator.py", line 138, in run 
 > check_genetic_models( 
 > File "/opt/conda/lib/python3.11/site-packages/genmod/annotate_models/genetic_models.py", line 195, in check_genetic_models 
 > elif check_compounds(variant_1, variant_2, family, intervals, phased): 
 > ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 
 > File "/opt/conda/lib/python3.11/site-packages/genmod/annotate_models/models/compound_model.py", line 86, in check_compounds 
 > variant_1_interval = intervals[individual_id].find_range( 
 > ~~~~~~~~~^^^^^^^^^^^^^^^ 
 > KeyError: 'HG002'

[dnil]

Hm, the current test cases for trios work (e.g. genmod annotate test_phased_compounds.vcf --annotate_regions | genmod models - -f recessive_trio.ped --phased) , so it's not really that it never works, rather that something goes wrong in the phase block calculation, possibly some chromosome with limited informative variants? I do think the approach with allowing precalculated PS tags is very prudent.

[dnil]

I can not reproduce this given the input data? The file given above appears to be without region annotation, so the command above is a bit incomplete, but I'm sure this is correct in the pipeline.

genmod annotate test_data_chunk.00002.bed_snv_genmod_annotate_annotate.vcf \
    --annotate_regions > test_data_chunk.00002.bed_snv_genmod_annotate_annotate.reg.vcf
genmod models --temp_dir . --phased --family_file FAM.ped \
    --processes 2 \
    --outfile test_data_chunk.00002.bed_snv_genmod_models_models.vcf \
    test_data_chunk.00002.bed_snv_genmod_annotate_annotate.reg.vcf

(genmod) ➜  phased_issue wc -l test_data_chunk.00002.bed_snv_genmod_models_models.vcf
     128 test_data_chunk.00002.bed_snv_genmod_models_models.vcf

Also, it seems to work just fine with the present VEP annotations, but then the --vep flag must ofcourse be used:

 genmod models --temp_dir . --vep --phased --family_file FAM.ped \
    --processes 2 \
    --outfile test_data_chunk.00002.bed_snv_genmod_models_models.vcf \
    test_data_chunk.00002.bed_snv_genmod_annotate_annotate.vcf

[fellen31]

I can not reproduce this given the input data? The file given above appears to be without region annotation, so the command above is a bit incomplete, but I'm sure this is correct in the pipeline.

How do you mean that test_data_chunk.00002.bed_snv_genmod_annotate_annotate.vcf.txt is without region annotation? It contains ##INFO=<ID=Annotation,Number=.,Type=String,Description="Annotates what feature(s) this variant belongs to."> and Annotation= for every variant?

[dnil]

Ah, some had! I just saw some didn't! 😊

genmod) ➜  phased_issue grep -v Annotation test_data_chunk.00002.bed_snv_genmod_annotate_annotate.vcf | grep -v \#
chr16	167947	chr16_167947_T_TC	T	TC	29	.	AF=1;AQ=29;FOUND_IN=deepvariant;cadd_raw=-2;cadd_phred=-2;gnomad_af=0.969714;CSQ=C|downstream_gene_variant|MODIFIER|HBM|ENSG00000206177|Transcript|ENST00000356815|protein_coding|||||||||||1183|1||HGNC|HGNC:4826|YES|1|P1|CCDS32347.1|ENSP00000349270|Q6B0K9.143|A0A1K0FU50.31|UPI000003730B|||||||||||||||||||||0.01,C|downstream_gene_variant|MODIFIER|HBM|ENSG00000206177|Transcript|ENST00000472539|protein_coding_CDS_not_defined|||||||||||1186|1||HGNC|HGNC:4826||5|||||||||||||||||||||||||||0.01,C|downstream_gene_variant|MODIFIER|HBM|ENSG00000206177|Transcript|ENST00000496585|protein_coding_CDS_not_defined|||||||||||1186|1||HGNC|HGNC:4826||2|||||||||||||||||||||||||||0.01;most_severe_consequence=HGNC:4826:C|downstream_gene_variant;most_severe_pli=0.01	GT:DP:AD:GQ:PL:RNC	1/1:15:5,9:8:29,7,0:..	1/1:9:1,7:21:29,21,0:..	1/1:9:2,7:21:24,23,0:..
chr16	168533	chr16_168533_T_C	T	C	33	.	AF=0.333333;AQ=33;FOUND_IN=deepvariant;cadd_raw=-2;cadd_phred=-2;gnomad_af=0.20605;CSQ=C|downstream_gene_variant|MODIFIER|HBM|ENSG00000206177|Transcript|ENST00000356815|protein_coding|||||||||||1769|1||HGNC|HGNC:4826|YES|1|P1|CCDS32347.1|ENSP00000349270|Q6B0K9.143|A0A1K0FU50.31|UPI000003730B|||||||||||||||||||||0.01,C|downstream_gene_variant|MODIFIER|HBM|ENSG00000206177|Transcript|ENST00000472539|protein_coding_CDS_not_defined|||||||||||1772|1||HGNC|HGNC:4826||5|||||||||||||||||||||||||||0.01,C|downstream_gene_variant|MODIFIER|HBM|ENSG00000206177|Transcript|ENST00000496585|protein_coding_CDS_not_defined|||||||||||1772|1||HGNC|HGNC:4826||2|||||||||||||||||||||||||||0.01;most_severe_consequence=HGNC:4826:C|downstream_gene_variant;most_severe_pli=0.01	GT:DP:AD:GQ:PL:RNC:PS	0|1:15:5,10:29:30,0,36:..:162329	0/0:9:9,0:27:0,27,269:..:.	1|0:10:6,4:33:33,0,41:..:160017
chr16	168736	chr16_168736_A_G	A	G	42	.	AF=0.666667;AQ=42;FOUND_IN=deepvariant;cadd_raw=-2;cadd_phred=-2;gnomad_af=0.772718;CSQ=G|downstream_gene_variant|MODIFIER|HBM|ENSG00000206177|Transcript|ENST00000356815|protein_coding|||||||||||1972|1||HGNC|HGNC:4826|YES|1|P1|CCDS32347.1|ENSP00000349270|Q6B0K9.143|A0A1K0FU50.31|UPI000003730B|||||||||||||||||||||0.01,G|downstream_gene_variant|MODIFIER|HBM|ENSG00000206177|Transcript|ENST00000472539|protein_coding_CDS_not_defined|||||||||||1975|1||HGNC|HGNC:4826||5|||||||||||||||||||||||||||0.01,G|downstream_gene_variant|MODIFIER|HBM|ENSG00000206177|Transcript|ENST00000496585|protein_coding_CDS_not_defined|||||||||||1975|1||HGNC|HGNC:4826||2|||||||||||||||||||||||||||0.01;most_severe_consequence=HGNC:4826:G|downstream_gene_variant;most_severe_pli=0.01	GT:DP:AD:GQ:PL:RNC:PS	0|1:15:6,9:37:37,0,46:..:162329	0|1:9:4,5:41:41,0,51:..:160017	1/1:9:0,9:38:42,39,0:..:.

[dnil]

After annotate regions (second pass of, then) these have got annotations.

(genmod) ➜  phased_issue grep -v Annotation test_data_chunk.00002.bed_snv_genmod_annotate_annotate.reg.vcf | grep -cv \#
0

And they do have VEP annotations, so they are fine with --vep.

[fellen31]

After annotate regions (second pass of, then) these have got annotations.

You need to annotate with genome build 38 for these files 😉

[dnil]

😄 Yes, I see. Ok, good, then we have the issue!