Logistic regression burden tests

Updates
3

Hi Enrique,

Regressions in Hail 0.2 are generic, so there are no longer special command for burden tests. Rather, to do logistic gene burden, annotate and group variants by gene to form a gene-by-sample matrix of scores, and then apply logistic regression per gene. Here’s an example that uses the total number minor alleles as a sample’s gene score (for a biallelic dataset with random phenotype and covariates):

import hail as hl

interval_ht = hl.import_locus_intervals('genes.37.interval_list')

mt = hl.import_vcf('profile.vcf.bgz')
mt = mt.annotate_cols(y = hl.rand_bool(0.5), x1 = hl.rand_norm(), x2 = hl.rand_norm())
mt = hl.variant_qc(mt)
mt = mt.annotate_rows(gene = interval_ht[mt.locus].target)
mt = mt.filter_rows(hl.is_defined(mt.gene))

gene_mt = mt.group_rows_by(mt.gene).aggregate(
    mac = hl.agg.sum(
        hl.cond(mt.variant_qc.AF[1] <= 0.5,
                mt.GT.n_alt_alleles(),
                2 - mt.GT.n_alt_alleles())))

gene_mt = hl.logistic_regression(y=gene_mt.y, x=gene_mt.mac, covariates=[1, gene_mt.x1, gene_mt.x2], test='wald')
gene_mt.logreg.show()

Here’s what’s shown:

+----------+--------------+-----------------------+---------------+----------------+
| gene     |  logreg.beta | logreg.standard_error | logreg.z_stat | logreg.p_value |
+----------+--------------+-----------------------+---------------+----------------+
| str      |      float64 |               float64 |       float64 |        float64 |
+----------+--------------+-----------------------+---------------+----------------+
| ATP6V1E1 |  1.69214e-03 |           1.28280e-03 |   1.31910e+00 |    1.87134e-01 |
| CCT8L2   |  4.03789e-02 |           3.46446e-02 |   1.16552e+00 |    2.43808e-01 |
| CECR1    |  1.29431e-03 |           1.82113e-03 |   7.10716e-01 |    4.77260e-01 |
| CECR2    |  4.90592e-04 |           3.88126e-04 |   1.26400e+00 |    2.06230e-01 |
| CECR5    | -2.90466e-04 |           3.41414e-03 |  -8.50774e-02 |    9.32200e-01 |
| CECR6    |  7.20520e-03 |           7.78411e-03 |   9.25630e-01 |    3.54638e-01 |
| GAB4     |  3.14315e-03 |           2.69332e-03 |   1.16701e+00 |    2.43205e-01 |
| IL17RA   |  3.40054e-03 |           4.98370e-03 |   6.82333e-01 |    4.95028e-01 |
| OR11H1   | -2.85557e-02 |           5.40733e-02 |  -5.28093e-01 |    5.97435e-01 |
| POTEH    | -6.64950e-03 |           3.76138e-03 |  -1.76784e+00 |    7.70884e-02 |
+----------+--------------+-----------------------+---------------+----------------+

+-------------------------+----------------------+---------------------+
| logreg.fit.n_iterations | logreg.fit.converged | logreg.fit.exploded |
+-------------------------+----------------------+---------------------+
|                   int32 | bool                 | bool                |
+-------------------------+----------------------+---------------------+
|                       3 | true                 | false               |
|                       3 | true                 | false               |
|                       3 | true                 | false               |
|                       3 | true                 | false               |
|                       2 | true                 | false               |
|                       3 | true                 | false               |
|                       3 | true                 | false               |
|                       3 | true                 | false               |
|                       3 | true                 | false               |
|                       3 | true                 | false               |
+-------------------------+----------------------+---------------------+
showing top 10 rows