# Import GLnexus pVCF

**URL:** <https://discuss.hail.is/t/import-glnexus-pvcf/2306>\
**Category:** Science\
**Created:** [October 25, 2021, 10:45pm UTC](https://discuss.hail.is/t/import-glnexus-pvcf/2306 "2021-10-25T22:45:02Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![poyingfu](https://avatars.discourse-cdn.com/v4/letter/p/2acd7d/32.png) [@poyingfu](https://discuss.hail.is/u/poyingfu)\
**Post date:** [October 25, 2021, 10:45pm UTC](https://discuss.hail.is/t/import-glnexus-pvcf/2306/1 "2021-10-25T22:45:02Z")

</div>

Hi Hail team,

I’m using GLnexus to joint-call my gvcfs, and I would like to import my pVCF into Hail to do further filtering. But some of GT like “./1”, “./0” will show as NA in Hail. Is there any possible way we can config to deal with this problem? Thanks!

Best,  
Po-Ying

---

<div class="post-metadata">

**Author:** ![danking](https://yyz2.discourse-cdn.com/flex036/user_avatar/discuss.hail.is/danking/32/43_2.png) [@danking](https://discuss.hail.is/u/danking)\
**Post date:** [October 26, 2021, 1:23pm UTC](https://discuss.hail.is/t/import-glnexus-pvcf/2306/2 "2021-10-26T13:23:49Z")

</div>

Hey @poyingfu, thanks for your question.

How would you like to interpret a `.`? Should `./1` be `0/1` or `1/1`?

---

<div class="post-metadata">

**Author:** ![kyle](https://yyz2.discourse-cdn.com/flex036/user_avatar/discuss.hail.is/kyle/32/11_2.png) [@kyle](https://discuss.hail.is/u/kyle)\
**Post date:** [October 27, 2021, 4:36am UTC](https://discuss.hail.is/t/import-glnexus-pvcf/2306/3 "2021-10-27T04:36:52Z")

</div>

This has to do with “MONOALLELIC” variants where GLnexus is trying to figure out what do with variants where a sample only has a copy number of 1 ([Reading GLnexus pVCFs · dnanexus-rnd/GLnexus Wiki · GitHub](https://github.com/dnanexus-rnd/GLnexus/wiki/Reading-GLnexus-pVCFs)).

I was working with the 2019 release of the SPARK data and ran my normal QC pipeline and ended up with a lot fewer indels than I expected. I looked into it and realized that all the genotypes at “MONOALLELIC” sites were being set to missing by Hail, so I had to figure out something to handle it.

I didn’t ultimately use the GLnexus version of the data (we ended up deciding to reprocess with GATK), but I found myself doing things like:

```auto
# This function looks for variants with "MONOALLELIC" in the filter field and:
# -If the genotype is missing, it makes it homref or het (homvar calls don't appear as missing)
# -Sets a value for the number of ref reads so that calls won't be filtered out later
def handle_monoallelics(mt):
    mt = mt.annotate_entries(GT = ( hl.case()
                               .when((mt.filters == {"MONOALLELIC"}) & hl.is_missing(mt.GT) & hl.is_missing(mt.AD[0]) & (mt.AD[1] == 0), 
                                    hl.Call([0,0]))
                               .when((mt.filters == {"MONOALLELIC"}) & hl.is_missing(mt.GT) & hl.is_missing(mt.AD[0]) & (mt.AD[1] > 0), 
                                    hl.Call([0,1]))
                               .default(mt.GT) ),
                           
                             AD = ( hl.case()
                               .when((mt.filters == {"MONOALLELIC"}) & hl.is_missing(mt.AD[0]), 
                                    [mt.DP - mt.AD[1], mt.AD[1]] )
                               .default(mt.AD) )
                            )
    return (mt)

```

This recovered my indels but I realized later it wasn’t doing everything I needed for variant call rate calculations, so I also did:

```auto
# This is to fix the no-calls where the RNC is ["-."] (a sample's call opposite its monoallelic variant)
mt = mt.annotate_entries(GT = hl.cond( (mt.RNC.size() == 1) & (mt.RNC[0] == "-.") & (hl.is_missing(mt.GT)) &
                                       (mt.AD[0] > 0) & (mt.AD[1] == 0), 
                                        hl.Call([0,0]), 
                                        mt.GT ),
                         
                         AD = hl.cond( (mt.RNC.size() == 1) & (mt.RNC[0] == "-.") & (hl.is_missing(mt.GT)) &
                                       (mt.AD[0] > 0) & (mt.AD[1] == 0), 
                                       [mt.DP, 0], 
                                        mt.AD )
                        )

```

So I don’t know if the format of the data you’re working with is exactly the same as this, but it’s example of how I tried to handle the issue you’re talking about.

---

<div class="post-metadata">

**Author:** ![poyingfu](https://avatars.discourse-cdn.com/v4/letter/p/2acd7d/32.png) [@poyingfu](https://discuss.hail.is/u/poyingfu)\
**Post date:** [October 28, 2021, 7:34pm UTC](https://discuss.hail.is/t/import-glnexus-pvcf/2306/4 "2021-10-28T19:34:58Z")

</div>

Thanks Kyle! That’s what I want to do!!! The way you approach it is pretty amazing.

My data look like below:

 ![Screen Shot 2021-10-28 at 1.14.15 PM](https://canada1.discourse-cdn.com/flex036/uploads/hail/original/1X/1533a18954e8cbc0e6367776195994df3d0db84a.jpeg)

 ![Screen Shot 2021-10-28 at 1.48.05 PM](https://canada1.discourse-cdn.com/flex036/uploads/hail/original/1X/18d20bba8ead7d3cdbca54708867b41aaab4af39.jpeg)

There are some cases, such as `./.`, `./0` and `./1`, I need to handle.

I want to do gene-base burden test, so I need to calculate hetero-ALT AC and homo-ALT AC for each variants.  
I’m really appreciate your help! It’s definitely a way I can try. Thanks a lot!