# Running polygenic risk scoring with VariantDataset

**URL:** <https://discuss.hail.is/t/running-polygenic-risk-scoring-with-variantdataset/3384>\
**Category:** Science\
**Created:** [May 15, 2023, 6:56am UTC](https://discuss.hail.is/t/running-polygenic-risk-scoring-with-variantdataset/3384 "2023-05-15T06:56:34Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![karenfeng](https://yyz2.discourse-cdn.com/flex036/user_avatar/discuss.hail.is/karenfeng/32/922_2.png) [@karenfeng](https://discuss.hail.is/u/karenfeng)\
**Post date:** [May 15, 2023, 6:56am UTC](https://discuss.hail.is/t/running-polygenic-risk-scoring-with-variantdataset/3384/1 "2023-05-15T06:56:34Z")

</div>

Do you have advice for how to run a polygenic risk score with a VariantDataset? My current flow is taking a painful amount of time to run on the All Of Us short-read WGS dataset and I’m not sure if it’s optimal.

```python
vds = hl.vds.read_vds(vds_wgs_path)
scores = hl.read_table(score_file)
vds = hl.vds.filter_samples(vds, samples, keep=True)
vds = hl.vds.split_multi(vds)
vds = hl.vds.filter_variants(vds, scores)
mt = hl.vds.to_dense_mt(vds)
mt = mt.annotate_rows(**scores[mt.locus, mt.alleles])
mt = hl.variant_qc(mt)
mt = mt.annotate_rows(prior=2 * hl.if_else(mt.variant_qc.AF.length() == 1, 0, mt.variant_qc.AF[1]))
mt = mt.annotate_cols(
prs=hl.agg.sum(
    mt['effect_allele_weight'] * hl.coalesce(mt.GT.n_alt_alleles(), mt.prior)))
mt.cols().export(prs_file)

```

---

<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:** [May 15, 2023, 1:21pm UTC](https://discuss.hail.is/t/running-polygenic-risk-scoring-with-variantdataset/3384/2 "2023-05-15T13:21:48Z")

</div>

AoU makes available an already-dense MatrixTable with just common variants, right? Using that as the starting point rather than the VDS (+ to\_dense\_mt) is probably going to be a lot more efficient for queries that don’t need rare variation.

---

<div class="post-metadata">

**Author:** ![danking](https://yyz2.discourse-cdn.com/flex036/user_avatar/discuss.hail.is/danking/32/43_2.png) [@danking](https://discuss.hail.is/u/danking)\
**Post date:** [May 15, 2023, 1:59pm UTC](https://discuss.hail.is/t/running-polygenic-risk-scoring-with-variantdataset/3384/3 "2023-05-15T13:59:23Z")

</div>

In addition to agreeing with Tim’s comments I’ll add that Hail does not yet have columnar storage (it’s in the works) so even running this operation on a dense matrix table will likely be bottlenecked by Hail’s decoding speed of data coming off the wire. This particularly bad for complex formats with lots of fields.

Also writing to a Hail Table will always yield better performance and avoid a concatenation step.

---

<div class="post-metadata">

**Author:** ![karenfeng](https://yyz2.discourse-cdn.com/flex036/user_avatar/discuss.hail.is/karenfeng/32/922_2.png) [@karenfeng](https://discuss.hail.is/u/karenfeng)\
**Post date:** [May 15, 2023, 4:43pm UTC](https://discuss.hail.is/t/running-polygenic-risk-scoring-with-variantdataset/3384/4 "2023-05-15T16:43:01Z")

</div>

Unfortunately, it seems like many of the variants are rare by AoU standards, so as a starting point I believe have to use the VDS. For this query, we’re using approximately 1M variants, which were common across a different population. I assume that the number of variants, along with the way they are semi-uniformly distributed, makes interval filtering a no-go.

I’m planning to pass these PRS to R for statistical analysis, so unfortunately writing to a text file may be the only way. Would it be cheaper to write it to a Hail table in between as a checkpointing step?

I noticed that `filter_samples` actually triggers a [row-wise filtering operation](https://github.com/hail-is/hail/blob/277ee08fc9cf6d8f972f66887f4ec4a430fae687/hail/python/hail/vds/methods.py#L393). Does this impede any interval filtering pushdown, so we should filter by interval first (if possible)? I was also surprised to see that interval filtering is expensive compared to MatrixTables, unless you specifically trigger a chromosome filter that allows you to use `unchecked_filter_both`.

This may be wishful thinking, but is it possible to just use the `variant_data` member of the VDS and get back into the MT regime without triggering `to_dense_mt`? Also, when would you recommend using a dense versus sparse MT?