# Exporting Hail MT to VCF - Missing Genotypes

**URL:** https://discuss.hail.is/t/exporting-hail-mt-to-vcf-missing-genotypes/3823
**Category:** Hail Query & hailctl
**Created:** [May 6, 2024, 5:16pm UTC](https://discuss.hail.is/t/exporting-hail-mt-to-vcf-missing-genotypes/3823 "2024-05-06T17:16:25Z")
**Posts on this page:** 12
**Page:** 1

<div class="post-metadata">

### Author: ![jasmussen](https://avatars.discourse-cdn.com/v4/letter/j/ad7895/32.png) [@jasmussen](https://discuss.hail.is/u/jasmussen)
#### Post date: [May 6, 2024, 5:16pm UTC](https://discuss.hail.is/t/exporting-hail-mt-to-vcf-missing-genotypes/3823/1 "2024-05-06T17:16:25Z")

</div>

Hi –

I am working in the Hail exome matrix table in the All of Us biobank. Once I subset the MT for my cohort of interest, I have been exporting to VCF. This is for a couple of reasons, but mainly b/c the algorithms used by my group read the variants from a VCF. In doing this, I’ve noticed that I lose about 30K variants per sample b/w the MT and VCF. Is this a known issue, or is it that the reference genotypes in the MT (0/0) are exported as missing in the VCF (./.)? If so, is there a way to bypass this issues, such that I get the reference calls in the VCF?

Thanks!

---

<div class="post-metadata">

### Author: ![patrick-schultz](https://yyz2.discourse-cdn.com/flex036/user_avatar/discuss.hail.is/patrick-schultz/32/265_2.png) [@patrick-schultz](https://discuss.hail.is/u/patrick-schultz)
#### Post date: [May 6, 2024, 7:14pm UTC](https://discuss.hail.is/t/exporting-hail-mt-to-vcf-missing-genotypes/3823/2 "2024-05-06T19:14:27Z")

</div>

Hi @jasmussen!  
I’m not sure what’s going on here. Could you share the code you’re running? If it’s possible to condense to a small example that shows the issue, that would make it easiest to identify the issue, but anything is helpful.

---

<div class="post-metadata">

### Author: ![jasmussen](https://avatars.discourse-cdn.com/v4/letter/j/ad7895/32.png) [@jasmussen](https://discuss.hail.is/u/jasmussen)
#### Post date: [May 6, 2024, 7:44pm UTC](https://discuss.hail.is/t/exporting-hail-mt-to-vcf-missing-genotypes/3823/3 "2024-05-06T19:44:20Z")

</div>

Hi @patrick-schultz!

Thank you for the response. Here is my basic approach. It is super simple -

First, I have a list of samples I want to include in the MT, so I filter columns with those IDs:  
mt = mt.filter\_cols(hl.literal(patient\_ids).contains(mt.s))

Then, I remove variants w/o genotypes in my set of samples:  
mt = hl.variant\_qc(mt)  
mt = mt.filter\_rows(mt.variant\_qc.AC[1] \> 0, keep = True)

Then, I export to VCF:  
output\_prefix = my\_bucket + ‘/outFile’  
hl.export\_vcf(mt, f’{output\_prefix}.vcf.bgz’)

That is it. When I run bcftools stats on the outFile.vcf.bgz I get about 30K less variants then what I started with in the Hail MT. I got the Hail MT per sample SNP counts by the following command:

mt = hl.sample\_qc(mt)  
ol\_df = mt.cols().to\_pandas()

Please let me know if you need anything else!

---

<div class="post-metadata">

### Author: ![patrick-schultz](https://yyz2.discourse-cdn.com/flex036/user_avatar/discuss.hail.is/patrick-schultz/32/265_2.png) [@patrick-schultz](https://discuss.hail.is/u/patrick-schultz)
#### Post date: [May 6, 2024, 8:03pm UTC](https://discuss.hail.is/t/exporting-hail-mt-to-vcf-missing-genotypes/3823/4 "2024-05-06T20:03:45Z")

</div>

> I got the Hail MT per sample SNP counts by the following command:

Was that run on the mt before or after filtering out variants with AC = 0?

I don’t know how bcftools computes stats, but if you reimport the vcf to a matrix table and compute the counts, do you still see a discrepency with before the export?

---

<div class="post-metadata">

### Author: ![jasmussen](https://avatars.discourse-cdn.com/v4/letter/j/ad7895/32.png) [@jasmussen](https://discuss.hail.is/u/jasmussen)
#### Post date: [May 7, 2024, 2:56am UTC](https://discuss.hail.is/t/exporting-hail-mt-to-vcf-missing-genotypes/3823/5 "2024-05-07T02:56:58Z")

</div>

I got the per sample SNP counts from the MT after filtering for AC = 0 variants.

Yes, I just tested it and after I re-import the VCF as a Hail MT the samples are no longer “missing” the 30K variants.

I thought ‘sample\_qc.n\_snp’ counted the total number of alternative alleles in each sample. I am very confused b/c bcftools stats should also count the total number of alternative alleles per sample (nNonRefHom + nHets). I’m not sure how to explain the difference.

---

<div class="post-metadata">

### Author: ![patrick-schultz](https://yyz2.discourse-cdn.com/flex036/user_avatar/discuss.hail.is/patrick-schultz/32/265_2.png) [@patrick-schultz](https://discuss.hail.is/u/patrick-schultz)
#### Post date: [May 7, 2024, 11:54am UTC](https://discuss.hail.is/t/exporting-hail-mt-to-vcf-missing-genotypes/3823/6 "2024-05-07T11:54:38Z")

</div>

Okay great, so the vcf isn’t missing any data, there’s just a discrepency in how the stats are being computed.

`sample_qc.n_snp` counts the number of alternate alleles _in snp sites_. I’m guessing the “missing” 30k variants are not snps. Try this to compute the count of all alternate alleles:

```python
mt.select_cols(n_alt_alleles=hl.agg.sum(mt.GT.n_alt_alleles)).cols().to_pandas()

```

---

<div class="post-metadata">

### Author: ![jasmussen](https://avatars.discourse-cdn.com/v4/letter/j/ad7895/32.png) [@jasmussen](https://discuss.hail.is/u/jasmussen)
#### Post date: [May 7, 2024, 5:14pm UTC](https://discuss.hail.is/t/exporting-hail-mt-to-vcf-missing-genotypes/3823/7 "2024-05-07T17:14:41Z")

</div>

Thank you for the feedback! That would make sense to me.

I was trying to run the command you provided and got an error. It ran after I include some additional parentheses:  
mt.select\_cols(n\_alt\_alleles=hl.agg.sum(mt.GT.n\_alt\_alleles())).cols().to\_pandas()

The output is a df with the sample ID and the sum of the n\_alt\_alleles, which is approx 122K per sample. That is 20K more than the sample\_qc.n\_snp. Is that what you were expecting?

---

<div class="post-metadata">

### Author: ![patrick-schultz](https://yyz2.discourse-cdn.com/flex036/user_avatar/discuss.hail.is/patrick-schultz/32/265_2.png) [@patrick-schultz](https://discuss.hail.is/u/patrick-schultz)
#### Post date: [May 7, 2024, 5:20pm UTC](https://discuss.hail.is/t/exporting-hail-mt-to-vcf-missing-genotypes/3823/8 "2024-05-07T17:20:03Z")

</div>

I was hoping it would arrive at the same counts that you’re getting from bcftools, to be sure there’s no discrepency in the exported vcf. It sounds like they still don’t quite agree?

---

<div class="post-metadata">

### Author: ![jasmussen](https://avatars.discourse-cdn.com/v4/letter/j/ad7895/32.png) [@jasmussen](https://discuss.hail.is/u/jasmussen)
#### Post date: [May 7, 2024, 5:21pm UTC](https://discuss.hail.is/t/exporting-hail-mt-to-vcf-missing-genotypes/3823/9 "2024-05-07T17:21:54Z")

</div>

Correct, they still don’t agree!

VCF bcftools stats = ~70K SNPs per sample  
Hail MT sample\_qc stats = ~100K SNPs per sample  
Hail MT n\_alt\_alleles = ~122K SNPs per sample

---

<div class="post-metadata">

### Author: ![patrick-schultz](https://yyz2.discourse-cdn.com/flex036/user_avatar/discuss.hail.is/patrick-schultz/32/265_2.png) [@patrick-schultz](https://discuss.hail.is/u/patrick-schultz)
#### Post date: [May 7, 2024, 5:41pm UTC](https://discuss.hail.is/t/exporting-hail-mt-to-vcf-missing-genotypes/3823/10 "2024-05-07T17:41:44Z")

</div>

> [@jasmussen](#):
>
> bcftools stats should also count the total number of alternative alleles per sample (nNonRefHom + nHets)

Oh, I missed this! The code I shared counted the number of alternate alleles, meaning a hom non-ref genotype adds 2 to the sum. This will count the number of sites with at least one alt allele:

```python
mt.select_cols(n_non_hom_ref=hl.agg.count(mt.GT.is_non_ref())).cols().to_pandas()

```

Hopefully that agrees with bcftools!

---

<div class="post-metadata">

### Author: ![jasmussen](https://avatars.discourse-cdn.com/v4/letter/j/ad7895/32.png) [@jasmussen](https://discuss.hail.is/u/jasmussen)
#### Post date: [May 8, 2024, 3:22pm UTC](https://discuss.hail.is/t/exporting-hail-mt-to-vcf-missing-genotypes/3823/11 "2024-05-08T15:22:18Z")

</div>

Alright, I didn’t get a chance to start it until late yesterday. The numbers are closer, but still don’t match bcftools. The average per sample is 89K, so more than bcftools. However, I did go back to the bcftools stats and this time I did the adjustment from that side of things. I did this –  
2\*nNonRefHom + nHets  
With doubling the nNonRefHom count, I get the same number as the sample\_qc.n\_snp from Hail.

Here is a recap –  
bcftools stats (nNonRefHom + nHets) = 73.6K sites with alternative alleles per sample  
bcftools stats (2\*nNonRefHom + nHets) = 101.9K alternative alleles per sample  
Hail MT sample\_qc.n\_snp = 101.9K alternative alleles per sample  
Hail MT n\_alt\_alleles = 122K alternative alleles per sample  
Hail MT n\_non\_hom\_ref = 88.6 sites that are non-hom-ref per sample

At this point, I feel good about the sample\_qc.n\_snp versus bcftools. Thank you so much for helping to sort it out!

Regarding the difference with Hail MT n\_non\_hom\_ref and bcftools (nNonRefHom + nHets), is Hail counting SNVs and indels that occur at the same site? I know you mentioned this in one of your earlier posts and, I know in bcftools I’m only counting SNVs.

---

<div class="post-metadata">

### Author: ![patrick-schultz](https://yyz2.discourse-cdn.com/flex036/user_avatar/discuss.hail.is/patrick-schultz/32/265_2.png) [@patrick-schultz](https://discuss.hail.is/u/patrick-schultz)
#### Post date: [May 8, 2024, 3:29pm UTC](https://discuss.hail.is/t/exporting-hail-mt-to-vcf-missing-genotypes/3823/12 "2024-05-08T15:29:31Z")

</div>

Great, glad we could get to the bottom of it!

> [@jasmussen](#):
>
> Regarding the difference with Hail MT n\_non\_hom\_ref and bcftools (nNonRefHom + nHets), is Hail counting SNVs and indels that occur at the same site? I know you mentioned this in one of your earlier posts and, I know in bcftools I’m only counting SNVs.

Right, the code snippets I shared are counting all variants, not just SNVs. I didn’t realize bcftools only counts SNVs, but it all makes sense now!
