# How to generate Manhattan plot for gene(burden) hail matrix table after regression analysis?

**URL:** <https://discuss.hail.is/t/how-to-generate-manhattan-plot-for-gene-burden-hail-matrix-table-after-regression-analysis/2266>\
**Category:** Hail Query & hailctl\
**Created:** [September 26, 2021, 3:41pm UTC](https://discuss.hail.is/t/how-to-generate-manhattan-plot-for-gene-burden-hail-matrix-table-after-regression-analysis/2266 "2021-09-26T15:41:59Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![Abhishek](https://yyz2.discourse-cdn.com/flex036/user_avatar/discuss.hail.is/abhishek/32/570_2.png) [@Abhishek](https://discuss.hail.is/u/Abhishek)\
**Post date:** [September 26, 2021, 3:41pm UTC](https://discuss.hail.is/t/how-to-generate-manhattan-plot-for-gene-burden-hail-matrix-table-after-regression-analysis/2266/1 "2021-09-26T15:41:59Z")

</div>

I have a gene(burden) hail matrix table where I condense all variants into genes for easier computation.

```python
mt = mt.annotate_entries(
            mac_adj=hl.if_else(
                mt.variant_qc.AF[1] <= 0.5,
                mt.GT.n_alt_alleles(),
                2 - mt.GT.n_alt_alleles()
            )
        )
gene_mt = mt.group_rows_by(mt.symbol, mt.feature).aggregate(
            mac=hl.agg.sum(mt.mac_adj),
            mac_HC=hl.agg.sum(hl.if_else(mt.lof == 'HC', mt.mac_adj, 0))

```

However, `group_rows_by` doesn’t have option of annotating with `locus` attribute.  
Later on, when I do regression analysis on matrix table and get the result, I am unable to plot Manhattan plot as it requires locus attribute.

> AttributeError: Table instance has no field, method, or property ‘locus’  
> Hint: use ‘describe()’ to show the names of all data fields.

How do you get Manhattan plot for gene matrix table based analysis?

PS: QQ plot is working fine

---

<div class="post-metadata">

**Author:** ![kumarveerapen](https://yyz2.discourse-cdn.com/flex036/user_avatar/discuss.hail.is/kumarveerapen/32/355_2.png) [@kumarveerapen](https://discuss.hail.is/u/kumarveerapen)\
**Post date:** [September 26, 2021, 6:17pm UTC](https://discuss.hail.is/t/how-to-generate-manhattan-plot-for-gene-burden-hail-matrix-table-after-regression-analysis/2266/2 "2021-09-26T18:17:23Z")

</div>

Hi, @Abhishek . Usually with gene burden tests, we describe our results with a table than with a Manhattan plot. If you do want to create a Manhattan plot, you could perhaps use `chr#:start-position` for each gene for visualization.

---

<div class="post-metadata">

**Author:** ![Abhishek](https://yyz2.discourse-cdn.com/flex036/user_avatar/discuss.hail.is/abhishek/32/570_2.png) [@Abhishek](https://discuss.hail.is/u/Abhishek)\
**Post date:** [September 26, 2021, 6:21pm UTC](https://discuss.hail.is/t/how-to-generate-manhattan-plot-for-gene-burden-hail-matrix-table-after-regression-analysis/2266/3 "2021-09-26T18:21:43Z")

</div>

Hi Veerapen,

Yes, I am generating a table separately as well for the results.

Can you provide me the hail code to generate `chr#:start-position` for each gene as a column, so I can annotate it in my final hail table(results) or gene hail matrix table?

---

<div class="post-metadata">

**Author:** ![Abhishek](https://yyz2.discourse-cdn.com/flex036/user_avatar/discuss.hail.is/abhishek/32/570_2.png) [@Abhishek](https://discuss.hail.is/u/Abhishek)\
**Post date:** [September 27, 2021, 5:28am UTC](https://discuss.hail.is/t/how-to-generate-manhattan-plot-for-gene-burden-hail-matrix-table-after-regression-analysis/2266/4 "2021-09-27T05:28:42Z")

</div>

I came up with this(which does not seem to be the optimized version but it works)

```auto
var_t = mt.rows()
var_t = var_t.group_by(vep_symbol=var_t.symbol, vep_feature=var_t.feature).aggregate(
            num_variants=hl.agg.count(),
            gene_contig=hl.agg.take(var_t.locus.contig, 1)[0],
            gene_pos=hl.agg.min(var_t.locus.position),
            gene_ref=hl.agg.take(var_t.locus.dtype.reference_genome.name, 1)[0]
            )

gene_mt = mt.group_rows_by(mt.symbol, mt.feature).aggregate(
            mac=hl.agg.sum(mt.mac_adj)
            )

gene_mt = gene_mt.annotate_rows(
            locus=hl.locus(var_t[gene_mt.row_key].gene_contig, var_t[gene_mt.row_key].gene_pos,
                           var_t[gene_mt.row_key].gene_ref.take(1)[0])
        )

```

---

<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, 2021, 12:30pm UTC](https://discuss.hail.is/t/how-to-generate-manhattan-plot-for-gene-burden-hail-matrix-table-after-regression-analysis/2266/5 "2021-09-27T12:30:59Z")

</div>

I think your last approach works, but I’ll say a bit more info about why a manhattan plot isn’t really appropriate here. We plot GWAS results using genomic coordinate as the x-axis and see “skyscrapers” in the results because linkage disequilibrium results in correlated genotypes and therefore correlated p-values in nearby sites. If you plot results per gene in the same way (using start locus or end locus or whatever), you won’t see any of the same structure because the results for nearby genes should not be correlated.

---

<div class="post-metadata">

**Author:** ![Abhishek](https://yyz2.discourse-cdn.com/flex036/user_avatar/discuss.hail.is/abhishek/32/570_2.png) [@Abhishek](https://discuss.hail.is/u/Abhishek)\
**Post date:** [September 28, 2021, 5:49am UTC](https://discuss.hail.is/t/how-to-generate-manhattan-plot-for-gene-burden-hail-matrix-table-after-regression-analysis/2266/6 "2021-09-28T05:49:04Z")

</div>

What will be a better plot for gene matrix table downstream analysis after regression?

---

<div class="post-metadata">

**Author:** ![kumarveerapen](https://yyz2.discourse-cdn.com/flex036/user_avatar/discuss.hail.is/kumarveerapen/32/355_2.png) [@kumarveerapen](https://discuss.hail.is/u/kumarveerapen)\
**Post date:** [September 28, 2021, 2:36pm UTC](https://discuss.hail.is/t/how-to-generate-manhattan-plot-for-gene-burden-hail-matrix-table-after-regression-analysis/2266/7 "2021-09-28T14:36:15Z")

</div>

Perhaps a volcano plot where the x-axis would be the effect size (log odds ratio) and y-axis being the -log10(p-value). Was inspired by a presentation by @Siwei this AM who used a similar visualization to present her gene-based test results today.

---

<div class="post-metadata">

**Author:** ![Abhishek](https://yyz2.discourse-cdn.com/flex036/user_avatar/discuss.hail.is/abhishek/32/570_2.png) [@Abhishek](https://discuss.hail.is/u/Abhishek)\
**Post date:** [September 28, 2021, 5:19pm UTC](https://discuss.hail.is/t/how-to-generate-manhattan-plot-for-gene-burden-hail-matrix-table-after-regression-analysis/2266/8 "2021-09-28T17:19:20Z")

</div>

It sounds like a good plot. Do you have a reference image or link to the presentation by @Siwei ?  
I would like to run through by my team.