# How to stat AC for each population

**URL:** <https://discuss.hail.is/t/how-to-stat-ac-for-each-population/662>\
**Category:** Help \[0.1\]\
**Created:** [September 25, 2018, 2:27am UTC](https://discuss.hail.is/t/how-to-stat-ac-for-each-population/662 "2018-09-25T02:27:09Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![whu\_wade](https://avatars.discourse-cdn.com/v4/letter/w/f19dbf/32.png) [@whu\_wade](https://discuss.hail.is/u/whu_wade)\
**Post date:** [September 25, 2018, 2:27am UTC](https://discuss.hail.is/t/how-to-stat-ac-for-each-population/662/1 "2018-09-25T02:27:09Z")

</div>

When I import a vcf file, and annotate each sample with population, age, gender and so on. Then how to stat each variant AC accross population, e.g. for variant1, AC of AFR, AMR and ASJ(The expected results showed below). I think it’s not efficient to filter each population and then calculate.  
AC  
AFR AMR ASJ …  
variant1 10,2 5,5 3,20  
variant2 2,32 1,5 5,10

---

<div class="post-metadata">

**Author:** ![tpoterba](https://yyz2.discourse-cdn.com/flex036/user_avatar/discuss.hail.is/tpoterba/32/61_2.png) [@tpoterba](https://discuss.hail.is/u/tpoterba)\
**Post date:** [September 25, 2018, 11:22am UTC](https://discuss.hail.is/t/how-to-stat-ac-for-each-population/662/2 "2018-09-25T11:22:38Z")

</div>

> I think it’s not efficient to filter each population and then calculate

This should be fine:

```python
mt = mt.annotate_rows(
    callstats_afr = hl.agg.call_stats(hl.agg.filter(mt.pop == 'AFR', mt.GT), mt.alleles),
    callstats_eur = hl.agg.call_stats(hl.agg.filter(mt.pop == 'EUR', mt.GT), mt.alleles),
...)

```

You can also use the experimental aggregator group\_by functionality:

```python
mt = mt.annotate_rows(call_stats = hl.agg.group_by(mt.pop, hl.agg.call_stats(mt.GT, mt.alleles)))

```

---

<div class="post-metadata">

**Author:** ![whu\_wade](https://avatars.discourse-cdn.com/v4/letter/w/f19dbf/32.png) [@whu\_wade](https://discuss.hail.is/u/whu_wade)\
**Post date:** [September 27, 2018, 6:16am UTC](https://discuss.hail.is/t/how-to-stat-ac-for-each-population/662/3 "2018-09-27T06:16:59Z")

</div>

Thanks. I have tried it on a test vcf, the remain question is that it cannot get/calculate the heterozygosis genotype counts when the allele number is over than 3

---

<div class="post-metadata">

**Author:** ![tpoterba](https://yyz2.discourse-cdn.com/flex036/user_avatar/discuss.hail.is/tpoterba/32/61_2.png) [@tpoterba](https://discuss.hail.is/u/tpoterba)\
**Post date:** [September 27, 2018, 1:48pm UTC](https://discuss.hail.is/t/how-to-stat-ac-for-each-population/662/4 "2018-09-27T13:48:27Z")

</div>

you can compute this separately:

```auto
mt = mt.annotate_rows(het_count = hl.agg.count_where(mt.GT.is_het()))

```

If you want to filter to AN \> 3:

```auto
mt = mt.filter_rows(mt.call_stats.AN > 3) # or what ever population

```

---

<div class="post-metadata">

**Author:** ![whu\_wade](https://avatars.discourse-cdn.com/v4/letter/w/f19dbf/32.png) [@whu\_wade](https://discuss.hail.is/u/whu_wade)\
**Post date:** [September 27, 2018, 3:09pm UTC](https://discuss.hail.is/t/how-to-stat-ac-for-each-population/662/5 "2018-09-27T15:09:47Z")

</div>

> [@tpoterba](#):
>
> hl.agg.group\_by(mt.pop

Maybe I haven’t express clear. For example there are four alleles in a site, ‘A’ for ref, and [‘G’,‘C’,‘T’] for alternative bases, the homo-genotype counts can be cal by call\_stats function, then how can I get the follow genotypes counts, AG, AC, AT …? Expect your reply, thanks

---

<div class="post-metadata">

**Author:** ![tpoterba](https://yyz2.discourse-cdn.com/flex036/user_avatar/discuss.hail.is/tpoterba/32/61_2.png) [@tpoterba](https://discuss.hail.is/u/tpoterba)\
**Post date:** [September 27, 2018, 4:02pm UTC](https://discuss.hail.is/t/how-to-stat-ac-for-each-population/662/6 "2018-09-27T16:02:23Z")

</div>

`annotate_rows(GC = hl.agg.counter(mt.GT))` should do it.

---

<div class="post-metadata">

**Author:** ![jjfarrell](https://yyz2.discourse-cdn.com/flex036/user_avatar/discuss.hail.is/jjfarrell/32/108_2.png) [@jjfarrell](https://discuss.hail.is/u/jjfarrell)\
**Post date:** [December 29, 2018, 5:23pm UTC](https://discuss.hail.is/t/how-to-stat-ac-for-each-population/662/7 "2018-12-29T17:23:39Z")

</div>

Similar to this, I tried call\_stats function for case and control stats but I am getting this assertion error. Any suggestions?

```
/restricted/projectnb/ukbiobank/ad/analysis/ad.v1/ad_parent_apoe_gt.py in <module>()
 73 mt = mt.filter_rows(mt.variant_qc.AF[0] <= (1.0-maf_cutoff))
 74 mt = mt.annotate_rows(
---> 75 case = hl.agg.call_stats(hl.agg.filter(mt.pheno.parent_AD_binary == 1, mt.GT), mt.alleles),
 76 control = hl.agg.call_stats(hl.agg.filter(mt.pheno.parent_AD_binary == 0, mt.GT), mt.alleles)
 77 )

<decorator-gen-995> in filter(condition, aggregation)

/restricted/projectnb/genpro/github/hail/hail/build/distributions/hail-python.zip/hail/typecheck/check.py in wrapper(__original_func, *args, **kwargs)
    557 @decorator
    558 def wrapper(__original_func, *args, **kwargs):
--> 559 args_, kwargs_ = check_all(__original_func, args, kwargs, checkers, is_method=is_method)
    560 return __original_func(*args_, **kwargs_)
    561

/restricted/projectnb/genpro/github/hail/hail/build/distributions/hail-python.zip/hail/typecheck/check.py in check_all(f, args, kwargs, checks, is_method)
    485 if i < len(args):
    486 arg = args[i]
--> 487 args_.append(checker.check(arg, name, arg_name))
    488 # passed as keyword
    489 else:

/restricted/projectnb/genpro/github/hail/hail/build/distributions/hail-python.zip/hail/expr/aggregators/aggregators.py in check(self, x, caller, param)
     24 def check(self, x, caller, param):
     25 x = self.coercer.check(x, caller, param)
---> 26 assert(len(x._ir.search(lambda node: isinstance(node, BaseApplyAggOp))) != 0)
     27 return x
     28

```

AssertionError:

---

<div class="post-metadata">

**Author:** ![tpoterba](https://yyz2.discourse-cdn.com/flex036/user_avatar/discuss.hail.is/tpoterba/32/61_2.png) [@tpoterba](https://discuss.hail.is/u/tpoterba)\
**Post date:** [December 29, 2018, 6:11pm UTC](https://discuss.hail.is/t/how-to-stat-ac-for-each-population/662/8 "2018-12-29T18:11:03Z")

</div>

ah, still need to fix that error message.

See here:

> [@\[BREAKING CHANGE\] redesign of aggregator interface](https://discuss.hail.is/t/breaking-change-redesign-of-aggregator-interface/701):
>
> the functions hl.agg.filter and hl.agg.explode have been changed. They now look like: hl.agg.filter(condition, hl.agg.collect(agg\_expr)) hl.agg.explode(lambda elt: hl.agg.collect(elt), agg\_array) They can be combined arbitrarily with each other and with hl.agg.group\_by, but must always an aggregation. For example: hl.agg.filter(t.idx \> 7, hl.agg.explode( lambda i: hl.agg.group\_by(t.group, hl.agg.collect(t.x)), …
