Logistic regression burden tests

Updates
1

We’ve added the method logreg_burden to provide a very general framework for aggregating (collapsing) genotypes and then applying logistic regression (Wald, LRT, Firth, or Score). This method returns both a key table of logistic regression results and a key table of collapsed sample scores for those samples with non-missing phenotype and covariates. The latter can then be further analyzed in PySpark or locally in tools like Python and R. Click the method link above for more details.

The performance is very similar to that of linreg_burden since there are typically far fewer keys for regression than variants for aggregation; it takes just a few minutes to run logistic gene burden starting from the annotated 2353 whole genomes in the 1000 Genomes Project.

2

Hi @jbloom,

Could you please point me out to a similar implementation of logreg_burden in Hail devel (0.2) version?

Thanks in advance!

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
3 Likes
4

Hi @jbloom

Thank for your quick reply. Could I aggregate by variant type (e.g. silent, damaging, etc…) and perform a logistic regression based on this same example? It looks like is possible to adapt it for performing a more global burden test using log_regression.

Thanks in advance

5

Absolutely, once you’ve annotated variants with their type, you can use group_rows_by and aggregate to create a type by sample matrix table (for example, the entries could count the number of variants of that type per sample) and then run logistic regression to test variant type.