# Filtering samples with extreme heterozygosity in hail?
**URL:** https://discuss.hail.is/t/filtering-samples-with-extreme-heterozygosity-in-hail/1277
**Category:** Hail Query & hailctl
**Created:** [February 6, 2020, 1:49pm UTC](https://discuss.hail.is/t/filtering-samples-with-extreme-heterozygosity-in-hail/1277 "2020-02-06T13:49:00Z")
**Posts on this page:** 8
**Page:** 1
### Author:  [@Raakhee](https://discuss.hail.is/u/Raakhee)
#### Post date: [February 6, 2020, 1:49pm UTC](https://discuss.hail.is/t/filtering-samples-with-extreme-heterozygosity-in-hail/1277/1 "2020-02-06T13:49:00Z")
Hi,
We are trying to use HAIL package to create a GWAS tutorial that is comparable to one we have made using R and Plink in hopes of being able to leverage Spark and possible improve the computational speed and/or memory consumption when conducting analyses (especially for larger datasets). Currently when attempting to perform quality control steps, we are able to get similar results when looking at variant call rate, minor allele frequency, and sample call rate outputs in both R and hail. When looking at heterozygosity to further assess sample quality and calculate a heterozygosity f statistic cut off to filter samples, we find sample qc statistics and available methods differ from R.
Our Current Approach in R:
The R snpStats package provides heterozygosity as an output which we plotted to determine what thresholds we wanted to use and manually calculated a heterozygosity f statistic/inbreeding cut off using minor allele frequency (for expected heterozygosity), and heterozygosity and N called for (observed heterozygosity). We then determined a f-statistic and calculated an f statistic cut off given the heterozygosity distribution to filter samples.
In Hail:
In hail we calculated mean heterozygosity by using the n\_het/n\_called columns output from the sample\_qc method.
What is the best approach for filtering samples with extreme heterozygosity in hail?
Currently we are attempting a similar approach to what we did in R and calculate a heterozygosity f statistic threshold by determining the max of the absolute value of the f statistic from samples within 2-2.5 standard deviations of the mean heterozygosity. To do this we plan to use the annotate\_cols and inbreeding aggregator to calculate the fstatistic and determine the thresholds using the n\_het/n\_called columns output from the sample qc method. Are there columns or stats included in the output from the sample qc method we should be using instead?
---
### Author:  [@tpoterba](https://discuss.hail.is/u/tpoterba)
#### Post date: [February 6, 2020, 2:02pm UTC](https://discuss.hail.is/t/filtering-samples-with-extreme-heterozygosity-in-hail/1277/2 "2020-02-06T14:02:23Z")
Can you provide the Hail Python script you’re running? I would expect this to work. I don’t have a sense of what snpStats is doing precisely, though.
---
### Author:  [@cseed](https://discuss.hail.is/u/cseed)
#### Post date: [February 6, 2020, 2:26pm UTC](https://discuss.hail.is/t/filtering-samples-with-extreme-heterozygosity-in-hail/1277/3 "2020-02-06T14:26:54Z")
An alternative to the heterozygosity f-statistic is to use the Hardy-Weinberg test p-value. This is computed by [hardy\_weinberg\_test](https://hail.is/docs/0.2/aggregators.html?highlight=hwe#hail.expr.aggregators.hardy_weinberg_test) aggregator and included in the statistics computed by the `variant_qc` method.
---
### Author:  [@Raakhee](https://discuss.hail.is/u/Raakhee)
#### Post date: [February 11, 2020, 4:26pm UTC](https://discuss.hail.is/t/filtering-samples-with-extreme-heterozygosity-in-hail/1277/4 "2020-02-11T16:26:53Z")
Thanks for the reply.
Below is the code for filtering samples with extreme heterozygosity in hail as well as R. Ultimately, the fstatistic value used as a cut off for filtering the samples with hail vs. R, but are similar. Given these results we have the following questions:
- Did we implement the filtering process correctly in hail? Is there any mistake we have made that could lead to the difference in the f-statistic cut off?
- If the hail code is correct should we consider the difference insignificant?
- If we were to investigate filtering samples using Hardy Weinberg as suggested by Cotton Seed below, how is the variant\_qc method used to filter samples? We do plan to filter variants using Hardy Weinberg as the next QC step, but were unsure how to apply this method when filtering samples.
> #calculate heterozygosity using the n\_het and n\_called fields from the sample\_qc method and fstatistic using the inbreeding method
>
> m1=m1.annotate\_cols(heterozygosity=((m1.sample\_qc.n\_het/m1.sample\_qc.n\_called)), inbreeding = hl.agg.inbreeding(m1.GT, m1.variant\_qc.AF[1]))
> #given the heterozygosity distribution, determine the min and the max heterozygosity values to be used when calculating a f-statistic cut off to filter samples
>
> het\_stats=m1.aggregate\_cols(hl.agg.stats(m1.heterozygosity))
>
> het\_low=(het\_stats.mean)-2.3\* (het\_stats.stdev)
>
> het\_high=(het\_stats.mean)+2.3 \*(het\_stats.stdev)
>
> het\_data= m1.filter\_cols((m1.heterozygosity \<= het\_high) & (m1.heterozygosity \>= het\_low))
> #calculate the fstatistic cut off
>
> het\_fstat\_cutoff = max(map(abs, het\_data.inbreeding.f\_stat.collect()))
> #filter samples filtered using fstatistic and sample call rate
>
> m1\_filtered\_fstat\_scallRate = m1.filter\_cols(((m1.inbreeding.f\_stat \<= abs(het\_fstat\_cutoff\_min)) &
>
> (m1.inbreeding.f\_stat \>= het\_fstat\_cutoff\_min)) &
>
> (m1.sample\_qc.call\_rate \> 0.95))
>
> m1\_filtered\_fstat\_scallRate.count\_cols()
The same approach followed in R is as follows:
> #heterozygosity is returned from the row.summary (snpsum.row) command from the snpStats package (we have compared distributions and they are they same)
>
> snpsum.row ← row.summary(genotype)
>
> callmatrix ← !is.na(genotype) _#get call matrix of complete calls_
> #calculate fstatististic (hetf) (using expected/observed heterozygosity)\*
>
> MAF ← snpsum.col$MAF
>
> hetExp ← callmatrix %_% (2_MAF\*(1-MAF))
>
> hetObs ← with(this\_snpsum.row,Heterozygosity\*(ncol(this\_genotype))\*Call.rate)
>
> snpsum.row$hetF ← 1-(hetObs/hetExp)
> #given the heterozygosity distribution, determine the min and the max heterozygosity values to be used when calculating a f-statistic cut off to filter samples
>
> snpsum.row\_low ← min(snpsum.row$Heterozygosity[snpsum.row$Heterozygosity \>= het\_low])
>
> snpsum.row\_up ← max(snpsum.row$Heterozygosity[snpsum.row$Heterozygosity \<= het\_up])
>
> snpsum.row\_up\_low ← snpsum.row[snpsum.row$Heterozygosity == snpsum.row\_low | snpsum.row$Heterozygosity == snpsum.row\_up,]
> #set heterozygosity F Statistic cut off
>
> hetcutoff ← max(abs(snpsum.row\_up\_low$hetF))
The outputs from both sets of code above are as follows
| **Calculation** | **Hail** | **R** |
| --- | --- | --- |
| Min heterozygosity value within in 2.3 sdev of mean | 0.31 | 0.31 |
| Max heterozygosity value within in 2.3 sdev of mean | 0.34 | 0.33 |
| F statistic cut off | 0.029 | 0.028 |
| Samples filtered | 33 (1368 remain) | 37 (1364 remain) |
When attempting to filter samples using Hardy Weinberg as Cotton Seed suggested, we tried the following:
> m1\_filtered\_hwe = m1.filter\_rows(m1.variant\_qc.p\_value\_hwe \> 10e-7)
> m1\_filtered\_hwe.count\_rows() #variants are filtered, yet the sample counts remain the same.
This makes sense given the method, but how should we be using it to filter samples?
---
### Author:  [@Raakhee](https://discuss.hail.is/u/Raakhee)
#### Post date: [February 12, 2020, 7:52pm UTC](https://discuss.hail.is/t/filtering-samples-with-extreme-heterozygosity-in-hail/1277/5 "2020-02-12T19:52:30Z")
Thanks for the reply.
I have mentioned the code in Hail and R. We would also want to understand few questions based on the code.
Below is the code for filtering samples with extreme heterozygosity in hail as well as R. Ultimately, the fstatistic value used as a cut off for filtering the samples with hail vs. R, but are similar. Given these results we have the following questions:
- Did we implement the filtering process correctly in hail? Is there any mistake we have made that could lead to the difference in the f-statistic cut off?
- If the hail code is correct should we consider the difference insignificant?
- If we were to investigate filtering samples using Hardy Weinberg as suggested by Cotton Seed below, how is the variant\_qc method used to filter samples? We do plan to filter variants using Hardy Weinberg as the next QC step, but were unsure how to apply this method when filtering samples.
> #calculate heterozygosity using the n\_het and n\_called fields from the sample\_qc method and fstatistic using the inbreeding method
>
> m1=m1.annotate\_cols(heterozygosity=((m1.sample\_qc.n\_het/m1.sample\_qc.n\_called)), inbreeding = hl.agg.inbreeding(m1.GT, m1.variant\_qc.AF[1]))
> #given the heterozygosity distribution, determine the min and the max heterozygosity values to be used when calculating a f-statistic cut off to filter samples
>
> het\_stats=m1.aggregate\_cols(hl.agg.stats(m1.heterozygosity))
>
> het\_low=(het\_stats.mean)-2.3\* (het\_stats.stdev)
>
> het\_high=(het\_stats.mean)+2.3 \*(het\_stats.stdev)
>
> het\_data= m1.filter\_cols((m1.heterozygosity \<= het\_high) & (m1.heterozygosity \>= het\_low))
> #calculate the fstatistic cut off
>
> het\_fstat\_cutoff = max(map(abs, het\_data.inbreeding.f\_stat.collect()))
> #filter samples filtered using fstatistic and sample call rate
>
> m1\_filtered\_fstat\_scallRate = m1.filter\_cols(((m1.inbreeding.f\_stat \<= abs(het\_fstat\_cutoff\_min)) &
>
> (m1.inbreeding.f\_stat \>= het\_fstat\_cutoff\_min)) &
>
> (m1.sample\_qc.call\_rate \> 0.95))
>
> m1\_filtered\_fstat\_scallRate.count\_cols()
The same approach followed in R is as follows:
> #heterozygosity is returned from the row.summary (snpsum.row) command from the snpStats package (we have compared distributions and they are they same)
>
> snpsum.row ← row.summary(genotype)
>
> callmatrix ← !is.na(genotype) _#get call matrix of complete calls_
> #calculate fstatististic (hetf) (using expected/observed heterozygosity)\*
>
> MAF ← snpsum.col$MAF
>
> hetExp ← callmatrix % _% (2_ MAF\*(1-MAF))
>
> hetObs ← with(this\_snpsum.row,Heterozygosity\*(ncol(this\_genotype))\*Call.rate)
>
> snpsum.row$hetF ← 1-(hetObs/hetExp)
> #given the heterozygosity distribution, determine the min and the max heterozygosity values to be used when calculating a f-statistic cut off to filter samples
>
> snpsum.row\_low ← min(snpsum.row$Heterozygosity[snpsum.row$Heterozygosity \>= het\_low])
>
> snpsum.row\_up ← max(snpsum.row$Heterozygosity[snpsum.row$Heterozygosity \<= het\_up])
>
> snpsum.row\_up\_low ← snpsum.row[snpsum.row$Heterozygosity == snpsum.row\_low | snpsum.row$Heterozygosity == snpsum.row\_up,]
> #set heterozygosity F Statistic cut off
>
> hetcutoff ← max(abs(snpsum.row\_up\_low$hetF))
The outputs from both sets of code above are as follows
| **Calculation** | **Hail** | **R** |
| --- | --- | --- |
| Min heterozygosity value within in 2.3 sdev of mean | 0.31 | 0.31 |
| Max heterozygosity value within in 2.3 sdev of mean | 0.34 | 0.33 |
| F statistic cut off | 0.029 | 0.028 |
| Samples filtered | 33 (1368 remain) | 37 (1364 remain) |
When attempting to filter samples using Hardy Weinberg as Cotton Seed suggested, we tried the following:
> m1\_filtered\_hwe = m1.filter\_rows(m1.variant\_qc.p\_value\_hwe \> 10e-7)
> m1\_filtered\_hwe.count\_rows() #variants are filtered, yet the sample counts remain the same.
This makes sense given the method, but how should we be using it to filter samples?
---
### Author:  [@kumarveerapen](https://discuss.hail.is/u/kumarveerapen)
#### Post date: [February 13, 2020, 6:46pm UTC](https://discuss.hail.is/t/filtering-samples-with-extreme-heterozygosity-in-hail/1277/6 "2020-02-13T18:46:35Z")
What I think is going on is the way that the sample missingness is handled in Hail vs R.
In R, you are explicitly removing sample rows using `!is.na`
As for Hail, we are handling missingness much like python. When one of the column items is missing, the resulting computation will result in missing. So, in Hail, you are not filtering out as many samples as you would in R due to missingness and that is likely the reason that you are having differing outputs.
---
### Author:  [@Raakhee](https://discuss.hail.is/u/Raakhee)
#### Post date: [February 13, 2020, 9:25pm UTC](https://discuss.hail.is/t/filtering-samples-with-extreme-heterozygosity-in-hail/1277/7 "2020-02-13T21:25:18Z")
Thanks Kumar. That makes sense. But also had one more question regarding using the variant\_qc.p\_value\_hwe present as one of the fields in the matrix table. How should we filter the samples on this field for the heterozygosity.
---
### Author:  [@tpoterba](https://discuss.hail.is/u/tpoterba)
#### Post date: [February 19, 2020, 1:37pm UTC](https://discuss.hail.is/t/filtering-samples-with-extreme-heterozygosity-in-hail/1277/8 "2020-02-19T13:37:22Z")
this is a variant field, not a sample field. Could you describe in words what you want to do to filter samples using HWE?