I’m trying to annotate a vcf with another vcf (gnomAD annotations), and the steps I follow are:
annotation_vds = hc.import_vcf(“annotations_vcf.vcf”).split_multi()
sample_vds = hc.import_vcf(“sample_vcf.vcf”).split_multi()
annotated_vds = sample_vds.annotate_variants_vds(annotation_vds," va.gnomAD_Ex_AF = vds.info.gnomAD_Ex_AF[vds.aIndex-1]")
Where vds.info.gnomAD_Ex_AF is an array with as many positions as alleles.
But when the annotation file (with gnomAD annotations) doesn’t have multiallelic variants, the ‘aIndex’ field doesn’t exist in the ‘va’ annotations. However, if I run ‘annotation_vds.was_split()’ it returns true.
So I can’t find a way to annotate the sample vcf if the annotation vcf doesn’t have multiallelic variants and the info field does. I’ve tried the following (without success):
annotated_vds = sample_vds.annotate_variants_vds(annotation_vds,“va.test = if(! isMissing(va.aIndex)) va.info.AC[va.aIndex-1] else va.info.AC”)
annotated_vds = sample_vds.annotate_variants_vds(annotation_vds,“va.test = orElse(va.info.AC[va.aIndex-1],va.info.AC)”)
Is this normal behaviour and there’s a way to annotate without knowing whether the vcf has multiallelic variants? Or is it not the expected behaviour?
Thanks in advance,