# PLINK glm vs linreg3 **URL:** \ **Category:** Hail Query & hailctl\ **Created:** [June 6, 2020, 8:38pm UTC](https://discuss.hail.is/t/plink-glm-vs-linreg3/1441 "2020-06-06T20:38:16Z")\ **Posts on this page:** 18\ **Page:** 1
**Author:** ![jimhavrilla](https://yyz2.discourse-cdn.com/flex036/user_avatar/discuss.hail.is/jimhavrilla/32/400_2.png) [@jimhavrilla](https://discuss.hail.is/u/jimhavrilla)\ **Post date:** [June 6, 2020, 8:38pm UTC](https://discuss.hail.is/t/plink-glm-vs-linreg3/1441/1 "2020-06-06T20:38:16Z")
Hey, curious about trying linreg3 for GWAS on UK Biobank and wondering how it compares to PLINK 2.0’s glm with logistic and fallback Firth regression? ---
**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:** [June 6, 2020, 11:38pm UTC](https://discuss.hail.is/t/plink-glm-vs-linreg3/1441/2 "2020-06-06T23:38:42Z")
linreg3 was a Hail 0.1 function that ran regressions on blocks of variants to use BLAS3 rather than BLAS2 routines to fit the linear model, which improved performance significantly. Hail 0.2 (which should definitely be used instead of Hail 0.1, and has seen ~8 person-years of work since 0.1 was deprecated) includes this functionality in the [`linear_regression_rows`](https://hail.is/docs/0.2/methods/stats.html#hail.methods.linear_regression_rows) function when multiple `y` values are passed. However, I should note that although 0.2 is superior to 0.1 on almost every performance benchmark, `import_bgen` / linear regression is one of the only benchmarks where 0.2 is **slightly slower than 0.1**. This is due to the presence of an extra copy into Hail’s off-heap data representation, and the fact that our current infrastructure is incompatible with a special-cased BGEN=\>linear regression optimization we used in 0.1. I expect that in the next few months, we’ll be able to reclaim 0.1’s performance, or improve it by using a faster on-disk format (zlib decompression is super slow!). In order to write something similar to PLINK’s --glm, you’d need to explicitly use linear or logistic regression depending on the type of phenotype (quantitative or binary). We also don’t have stuff to do firth fallback right now, though I think that’s not a super hard feature to add. ---
**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:** [June 6, 2020, 11:44pm UTC](https://discuss.hail.is/t/plink-glm-vs-linreg3/1441/3 "2020-06-06T23:44:20Z")
Also, to answer what may be your original question, I believe that 0.1’s linreg3 and 0.2’s linear\_regression\_rows are faster than PLINK for multiple-phenotype regression, due to the BLAS3 optimizations I mentioned above. Hail will probably be a little slower (in terms of CPU time, Hail’s scalability means you can use a cluster more easily to parallelize compute) for logistic regression, though may be faster if you’re running on multiple phenotypes at once, since I think PLINK will need to go to disk for every phenotype, while Hail can batch them. ---
**Author:** ![jimhavrilla](https://yyz2.discourse-cdn.com/flex036/user_avatar/discuss.hail.is/jimhavrilla/32/400_2.png) [@jimhavrilla](https://discuss.hail.is/u/jimhavrilla)\ **Post date:** [June 7, 2020, 12:02am UTC](https://discuss.hail.is/t/plink-glm-vs-linreg3/1441/4 "2020-06-07T00:02:43Z")
Thanks for the thorough response. Clarification question, what do you mean PLINK has to do I/O while Hail can batch phenotype? I can set up qsub on Sun Grid Engine to run one phenotype at a time no problem. My main concern is for large bgens ~150GB per chrom, I can’t reasonably merge them into one file and perform PCA and logistic regression. I understand from Ben Neale’s work using your software it is quite feasible to run logreg using all of the chromosomes at once but I’m a bit confused on how it would incorporate them all. I think PLINK would need to run them one by one because the BGENs are so large or else hit OOM errors. Maybe I’m misunderstanding something? Perhaps a throwaway example command might give some clarification. I’m probably just not getting it. Jim Havrilla ---
**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:** [June 7, 2020, 12:24am UTC](https://discuss.hail.is/t/plink-glm-vs-linreg3/1441/5 "2020-06-07T00:24:52Z")
> what do you mean PLINK has to do I/O while Hail can batch phenotype? I’m speaking to computational efficiency, which is related to cost when you’re paying marginally per CPU hour. It sounds like that’s not your core concern, though. Hail’s scalability means that you can write the same code, whether you’re running on a 1GB dataset or 1T. You can process everything at once without having to worry about qsub/SGE! Roughly the code would look like the below: ```python import hail as hl list_of_bgens = [f'/path/to/chr{chrom}.bgen' for chrom in range(1, 23)] mt = hl.import_bgen(list_of_bgens, entry_fields=('GT', 'GP'), sample_file='/path/to/sample/file') # do any necessary sample/variant filtering here # join phenotypes into matrix table pheno_table = hl.import_table('/path/to/pheno_file.tsv', ...) mt = mt.annotate_cols(phenos = pheno_table[mt.s]) gwas_results = hl.linear_regression_rows( y=[mt.pheno.pheno1, mt.pheno.pheno2, ...], x=mt.GP, covariates=[1, ...pcs...]) ``` This will scale to the full biobank. Nothing in Hail assumes that you can put an entire dataset in memory – methods are out-of-core, running row-by-row through the data. ---
**Author:** ![jimhavrilla](https://yyz2.discourse-cdn.com/flex036/user_avatar/discuss.hail.is/jimhavrilla/32/400_2.png) [@jimhavrilla](https://discuss.hail.is/u/jimhavrilla)\ **Post date:** [June 8, 2020, 6:46pm UTC](https://discuss.hail.is/t/plink-glm-vs-linreg3/1441/6 "2020-06-08T18:46:14Z")
Thanks, and this is probably a stupid question (and possibly a bit unrelated), but is there a way to pass multiple BGENs in for PCA and logistic/linear regression for running on the whole genome at once instead of doing QC and merging into one unified file as with PLINK (at least that’s how I understand it)? ---
**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:** [June 8, 2020, 7:51pm UTC](https://discuss.hail.is/t/plink-glm-vs-linreg3/1441/7 "2020-06-08T19:51:37Z")
Yes. `import_bgen` accepts a list of files. ---
**Author:** ![jimhavrilla](https://yyz2.discourse-cdn.com/flex036/user_avatar/discuss.hail.is/jimhavrilla/32/400_2.png) [@jimhavrilla](https://discuss.hail.is/u/jimhavrilla)\ **Post date:** [June 9, 2020, 2:30am UTC](https://discuss.hail.is/t/plink-glm-vs-linreg3/1441/8 "2020-06-09T02:30:21Z")
Thanks! Good to know I can use multiple files for PCA and linreg at once! Another probably similar stupid question…can I pass multiple gVCFs in for a SKAT test? Or do I have to joint call them/merge using GLnexus or GATK4? Because currently merging them is a big bottleneck for me. ---
**Author:** ![cseed](https://yyz2.discourse-cdn.com/flex036/user_avatar/discuss.hail.is/cseed/32/361_2.png) [@cseed](https://discuss.hail.is/u/cseed)\ **Post date:** [June 9, 2020, 2:56am UTC](https://discuss.hail.is/t/plink-glm-vs-linreg3/1441/9 "2020-06-09T02:56:01Z")
You’ll want to merge first before running SKAT. However, Hail has a joint caller: [https://hail.is/docs/0.2/experimental/vcf\_combiner.html?highlight=vcf\_combiner](https://hail.is/docs/0.2/experimental/vcf_combiner.html?highlight=vcf_combiner). The output differs from GATK4, etc., it is a “sparse MatrixTable” format that preserves REF blocks and is designed to support very large (1m+ sample) datasets. You’ll want to densify after joint calling to get an analog of the traditional joint-called VCF format, see the docs. ---
**Author:** ![jimhavrilla](https://yyz2.discourse-cdn.com/flex036/user_avatar/discuss.hail.is/jimhavrilla/32/400_2.png) [@jimhavrilla](https://discuss.hail.is/u/jimhavrilla)\ **Post date:** [June 9, 2020, 3:39am UTC](https://discuss.hail.is/t/plink-glm-vs-linreg3/1441/10 "2020-06-09T03:39:17Z")
Ok, so I should not use import\_gvcfs (though the run did finish) Another question that was a bit unclear to me from the documentation here ([https://hail.is/docs/0.2/getting\_started.html#installation](https://hail.is/docs/0.2/getting_started.html#installation)). If I am running this on the head node of an SGE based cluster (I would assume that’s a spark cluster) is the Spark jobs page that I can access (in my case) from port 4040 based solely on the head node, or does Spark intuitively use the machine resources on the other compute nodes? Normally I have to qsub on my cluster and I know that spark-qsub and spark-qrsh exist, and I can’t seem to see how it’s distributing the memory usage with qstat. I just want to be certain I’m not going to be crashing the head node or anything like that if I run a full, non-test-size job – for now I’m just looking at one gene across all the gVCFs and it seems to load successfully. I also want to say you guys have been insanely helpful and I really appreciate it. ---
**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:** [June 9, 2020, 1:49pm UTC](https://discuss.hail.is/t/plink-glm-vs-linreg3/1441/11 "2020-06-09T13:49:37Z")
Spark and SGE are distinct cluster management systems. If a `spark-qsub` binary exists, then I imagine your IT department has somehow built a Spark-on-SGE system. Usually there’s a binary called `spark-submit` which can accept a python file which starts a Spark job. Can you `spark-qsub /path/to/hail-script.py`? We also maintain tools for using managed Spark clusters in [Google Cloud](https://hail.is/docs/0.2/cloud/google_cloud.html). You can also run Hail on one computer with many cores. Hail will use all the cores available on that machine. ---
**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:** [June 9, 2020, 1:55pm UTC](https://discuss.hail.is/t/plink-glm-vs-linreg3/1441/12 "2020-06-09T13:55:12Z")
Just to be a bit more explicit, try this: ```bash cat >my-script.py < **Author:** ![jimhavrilla](https://yyz2.discourse-cdn.com/flex036/user_avatar/discuss.hail.is/jimhavrilla/32/400_2.png) [@jimhavrilla](https://discuss.hail.is/u/jimhavrilla)\ **Post date:** [June 9, 2020, 6:20pm UTC](https://discuss.hail.is/t/plink-glm-vs-linreg3/1441/13 "2020-06-09T18:20:20Z") Yeah we don’t have spark-qsub set up. I was more curious how I could use hail on the cluster to make best use of my resources. My guess is when I run it, it is just using the head node. I guess I’m asking should I be starting an interactive session or qsub job so that Hail can have access to more resources when I run it? Or am I limiting it or perhaps breaking it if I do so? ---
**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:** [June 9, 2020, 6:22pm UTC](https://discuss.hail.is/t/plink-glm-vs-linreg3/1441/14 "2020-06-09T18:22:19Z")
You can definitely run on a beefy interactive node or in a single `qsub` job that requests a node with lots of resources. This will be running Spark in local mode, which will allocate all cores it finds by default. You can set the cores used by using the following in your python script: ```auto import hail as hl hl.init(master='local[N]') ``` ---
**Author:** ![jimhavrilla](https://yyz2.discourse-cdn.com/flex036/user_avatar/discuss.hail.is/jimhavrilla/32/400_2.png) [@jimhavrilla](https://discuss.hail.is/u/jimhavrilla)\ **Post date:** [June 9, 2020, 6:26pm UTC](https://discuss.hail.is/t/plink-glm-vs-linreg3/1441/15 "2020-06-09T18:26:03Z")
Caveat being that I ran it on the head node and it is able to create matrix tables for a small gene in 49k VCFs no problem (still working on fixing vcf\_combiner though), but I haven’t tried the whole exome yet, so I don’t know how well that will work out. ---
**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:** [June 9, 2020, 6:30pm UTC](https://discuss.hail.is/t/plink-glm-vs-linreg3/1441/16 "2020-06-09T18:30:17Z")
> [@jimhavrilla](#): > > Yeah we don’t have spark-qsub set up. I was more curious how I could use hail on the cluster to make best use of my resources. My guess is when I run it, it is just using the head node. I guess I’m asking should I be starting an interactive session or qsub job so that Hail can have access to more resources when I run it? Or am I limiting it or perhaps breaking it if I do so? These questions are best answered by the IT department that supports and maintains your SGE cluster. I’ve never seen Spark running on SGE before, so I suspect whatever they did is custom to your institution. ---
**Author:** ![jimhavrilla](https://yyz2.discourse-cdn.com/flex036/user_avatar/discuss.hail.is/jimhavrilla/32/400_2.png) [@jimhavrilla](https://discuss.hail.is/u/jimhavrilla)\ **Post date:** [June 9, 2020, 6:36pm UTC](https://discuss.hail.is/t/plink-glm-vs-linreg3/1441/17 "2020-06-09T18:36:01Z")
Awesome, thanks! Once I find a quick workaround for merging those gVCF matrix tables into one, I’m in business. Jim Havrilla ---
**Author:** ![jimhavrilla](https://yyz2.discourse-cdn.com/flex036/user_avatar/discuss.hail.is/jimhavrilla/32/400_2.png) [@jimhavrilla](https://discuss.hail.is/u/jimhavrilla)\ **Post date:** [June 9, 2020, 6:41pm UTC](https://discuss.hail.is/t/plink-glm-vs-linreg3/1441/18 "2020-06-09T18:41:51Z")
Thanks I think Tim answered my question pretty well. I do want to thank you all again for being very polite and trying to help me out. It’s definitely appreciated. Jim Havrilla