# Fast linear regression for multiple phenotypes, eQTLs **URL:** \ **Category:** Updates\ **Created:** [April 16, 2017, 3:09pm UTC](https://discuss.hail.is/t/fast-linear-regression-for-multiple-phenotypes-eqtls/190 "2017-04-16T15:09:34Z")\ **Posts on this page:** 20\ **Page:** 1
**Author:** ![jbloom](https://yyz2.discourse-cdn.com/flex036/user_avatar/discuss.hail.is/jbloom/32/109_2.png) [@jbloom](https://discuss.hail.is/u/jbloom)\ **Post date:** [April 16, 2017, 3:09pm UTC](https://discuss.hail.is/t/fast-linear-regression-for-multiple-phenotypes-eqtls/190/1 "2017-04-16T15:09:34Z")
We’ve added a version of linear regression called [linreg\_multi\_pheno](https://hail.is/hail/hail.VariantDataset.html?highlight=linreg_multi_pheno#hail.VariantDataset.linreg_multi_pheno) that takes a list of phenotypes rather than a single phenotype, and is far more efficient than looping over [linreg](https://hail.is/hail/hail.VariantDataset.html?highlight=linreg#hail.VariantDataset.linreg) for scenarios where the set of samples is fixed over all regressions. In particular, this method only makes one pass over the data and takes greater advantage of BLAS optimization. We took the method for a one-billion-regression-[eQTL](https://en.wikipedia.org/wiki/Expression_quantitative_trait_loci)-spin on a Macbook Pro. Namely: 20k phenotypes (gene expression from [Geuvadis](http://www.geuvadis.org)), 50k common variants on 464 sample from the 1000 genomes project, with 15 covariates (intercept, sex, 3 PCs, 10 [Peer factors](http://www.sanger.ac.uk/science/tools/peer)). The billion multivariate regressions took about 150 seconds, whereas exporting the p-values to a 9GB compressed tsv.gz took 8 minutes more. Instead, pushing the p-values to a PySpark Dataframe and saving to Parquet format only took 1 minute more (although it was still 8GB). Of course, one might prefer to just continue analysis in PySpark and write a far smaller file of interesting results. And the same method may be run variant-parallel on thousands of cores in the [cloud](http://discuss.hail.is/t/using-hail-on-the-google-cloud-platform/80), making whole-genome cis and trans eQTL studies doable on the order of minutes. ---
**Author:** ![davidkelley](https://avatars.discourse-cdn.com/v4/letter/d/6bbea6/32.png) [@davidkelley](https://discuss.hail.is/u/davidkelley)\ **Post date:** [September 26, 2017, 6:03pm UTC](https://discuss.hail.is/t/fast-linear-regression-for-multiple-phenotypes-eqtls/190/2 "2017-09-26T18:03:37Z")
Hey Jon, This sounds really cool! If one wanted to use more variants, but focus on a cis-region around each gene, would that play nice with the existing Hail framework? I’ve run some initial experiments where I naively use filter\_intervals, but that seems to slow things down considerably relative to the performance you’re seeing here. Best, Dave ---
**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, 2017, 3:58pm UTC](https://discuss.hail.is/t/fast-linear-regression-for-multiple-phenotypes-eqtls/190/3 "2017-09-27T15:58:20Z")
@DylanVrana did something similar (eQTLs in 1MB window around each gene) – I’m not sure how exactly his stuff went, but he would be able to say more. ---
**Author:** ![DylanVrana](https://avatars.discourse-cdn.com/v4/letter/d/b782af/32.png) [@DylanVrana](https://discuss.hail.is/u/DylanVrana)\ **Post date:** [September 27, 2017, 6:08pm UTC](https://discuss.hail.is/t/fast-linear-regression-for-multiple-phenotypes-eqtls/190/4 "2017-09-27T18:08:44Z")
Hi @davidkelley, I have been using filter\_intervals for my cis-eQTL-style analysis (across GTEx tissues, so ~36 bn regressions). I’m not sure of the performance costs of filter\_intervals (@tpoterba probably knows better than I do- any thoughts about it?). My filter\_intervals based cis-eQTL analyses take between 30 minutes and 2 hours per tissue (~1.2 bn regressions, 100-200 samples) on a 6-worker Dataproc cluster, depending on the number of samples and how often the workers get preempted away. ---
**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, 2017, 7:08pm UTC](https://discuss.hail.is/t/fast-linear-regression-for-multiple-phenotypes-eqtls/190/5 "2017-09-27T19:08:37Z")
Can you share your code to do that? I kinda remember you doing some `union` stuff to parallelize among intervals. ---
**Author:** ![jerryliu2005](https://yyz2.discourse-cdn.com/flex036/user_avatar/discuss.hail.is/jerryliu2005/32/97_2.png) [@jerryliu2005](https://discuss.hail.is/u/jerryliu2005)\ **Post date:** [September 29, 2017, 1:04pm UTC](https://discuss.hail.is/t/fast-linear-regression-for-multiple-phenotypes-eqtls/190/6 "2017-09-29T13:04:59Z")
We really like the linreg\_multi\_pheno function. Is there an equivalent logreg function for multiple phenotypes with categorical values in the cooking or already done? Thanks! ---
**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:** [October 3, 2017, 5:14pm UTC](https://discuss.hail.is/t/fast-linear-regression-for-multiple-phenotypes-eqtls/190/7 "2017-10-03T17:14:20Z")
Sorry, missed this last week. We should definitely have this, and I think it’s something we could expect to be in 0.2. Part of the reason for linreg\_multi\_pheno is the speedup from vectorizing the matrix multiplies, though, and logistic regression won’t benefit from that boost. ---
**Author:** ![jerryliu2005](https://yyz2.discourse-cdn.com/flex036/user_avatar/discuss.hail.is/jerryliu2005/32/97_2.png) [@jerryliu2005](https://discuss.hail.is/u/jerryliu2005)\ **Post date:** [October 5, 2017, 1:05pm UTC](https://discuss.hail.is/t/fast-linear-regression-for-multiple-phenotypes-eqtls/190/8 "2017-10-05T13:05:17Z")
Thanks. When is 0.2 expected to be out? ---
**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:** [October 5, 2017, 2:20pm UTC](https://discuss.hail.is/t/fast-linear-regression-for-multiple-phenotypes-eqtls/190/9 "2017-10-05T14:20:13Z")
Can’t really say – there’s still a lot to do. My rough guess would be 3-6 months, but that was also my guess 3 months ago. ---
**Author:** ![Stephane\_Bourgeois](https://avatars.discourse-cdn.com/v4/letter/s/f04885/32.png) [@Stephane\_Bourgeois](https://discuss.hail.is/u/Stephane_Bourgeois)\ **Post date:** [October 6, 2017, 10:20pm UTC](https://discuss.hail.is/t/fast-linear-regression-for-multiple-phenotypes-eqtls/190/10 "2017-10-06T22:20:40Z")
Hi @DylanVrana, it would be really great if you could share how you feed a (large) list of intervals to filter\_intervals. Thanks, Steph ---
**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:** [October 7, 2017, 1:15am UTC](https://discuss.hail.is/t/fast-linear-regression-for-multiple-phenotypes-eqtls/190/11 "2017-10-07T01:15:54Z")
Dylan’s use case was a bit special – he was doing cis-eQTL analysis, so he wanted to run regressions for variants near each gene. The natural way to do this was looping through genes, and: 1. Annotation the expression for that gene 2. Filtering the dataset to intervals near that gene 3. Running linear regression and exporting results. Dylan and I worked together to formulate a better strategy for this problem. If you just want to filter to a dataset to a large set of intervals (say, filtering WGS data to exome capture regions), then you should use [`filter_variants_table`](https://hail.is/docs/stable/hail.VariantDataset.html#hail.VariantDataset.filter_variants_table). You can import an interval list into a key table with [`KeyTable.import_interval_list`](https://hail.is/docs/stable/hail.KeyTable.html#hail.KeyTable.import_interval_list) or [`KeyTable.import_bed`](https://hail.is/docs/stable/hail.KeyTable.html#hail.KeyTable.import_bed). ---
**Author:** ![Stephane\_Bourgeois](https://avatars.discourse-cdn.com/v4/letter/s/f04885/32.png) [@Stephane\_Bourgeois](https://discuss.hail.is/u/Stephane_Bourgeois)\ **Post date:** [October 7, 2017, 10:16am UTC](https://discuss.hail.is/t/fast-linear-regression-for-multiple-phenotypes-eqtls/190/12 "2017-10-07T10:16:05Z")
These functions will indeed be very useful when filtering exomes, thanks. What I’m also after is similar to Dylan’s situation, I’m trying to find out if I can analyse methylation data in Hail using univariate or multivariate linear mixed models (similar to what GEMMA does [http://www.xzlab.org/software.html](http://www.xzlab.org/software.html)). We’re doing cis and trans analyses, and the way GEMMA works, we have to create hundreds of thousands of files, one for each set of variants for which we want to perform the analysis. The way Hail works would alleviate this burden, and allow us to perform the analyses by filtering based on intervals. Do you think that’s possible? Is there a way to annotate g with methylation status? ---
**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:** [October 7, 2017, 2:57pm UTC](https://discuss.hail.is/t/fast-linear-regression-for-multiple-phenotypes-eqtls/190/13 "2017-10-07T14:57:14Z")
The genotype field in 0.1 is generally expected to be the GATK schema: GT, AD, DP, GQ, PL. We do have limited support for totally generic genotypes, but you can’t do much with them and everything runs much more slowly on them. Many methods are hard-coded against this schema, like linear regression (it only regresses against GT or dosage). We expect to have a completely generic schema in 0.2 with no performance cost, which will let you annotate the methylation status and use that in linear regression, for instance. ---
**Author:** ![dwerling](https://yyz2.discourse-cdn.com/flex036/user_avatar/discuss.hail.is/dwerling/32/138_2.png) [@dwerling](https://discuss.hail.is/u/dwerling)\ **Post date:** [April 11, 2018, 9:00pm UTC](https://discuss.hail.is/t/fast-linear-regression-for-multiple-phenotypes-eqtls/190/14 "2018-04-11T21:00:05Z")
Hi there - I am also working on a cis-eQTL analysis at the moment, but I keep getting a timeout error when I try to loop through all genes on a dataproc cluster. Did you and Dylan run into a similar issue, and if so, how did you solve it? Would appreciate any insight you have about running and optimizing an analysis like this. Thanks, Donna ---
**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:** [April 11, 2018, 9:29pm UTC](https://discuss.hail.is/t/fast-linear-regression-for-multiple-phenotypes-eqtls/190/15 "2018-04-11T21:29:38Z")
can you provide the pipeline and the error messages? ---
**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:** [April 11, 2018, 9:29pm UTC](https://discuss.hail.is/t/fast-linear-regression-for-multiple-phenotypes-eqtls/190/16 "2018-04-11T21:29:57Z")
This is Hail 0.2 right? ---
**Author:** ![dwerling](https://yyz2.discourse-cdn.com/flex036/user_avatar/discuss.hail.is/dwerling/32/138_2.png) [@dwerling](https://discuss.hail.is/u/dwerling)\ **Post date:** [April 11, 2018, 11:40pm UTC](https://discuss.hail.is/t/fast-linear-regression-for-multiple-phenotypes-eqtls/190/17 "2018-04-11T23:40:53Z")
Hi - no this is 0.1. I would switch, but we’re in the middle of trying to meet a tight deadline and I’m reluctant to try to learn the new version…though I could be convinced! Here’s the pipeline, as run on dataproc using Jupyter notebook (apologies in advance for length). The final block with the loop is where it falls apart: ```auto # Load Hail import hail from hail import * hc = HailContext() # Load packages import pandas as pd import numpy as np from math import log, isnan import os,sys,glob,gzip from pprint import pprint ``` ```auto # Load toy VDS vds = hc.read('gs://donnawerlingucsf/devEQTL/filterPart2_toy-180402.vds') ``` ```auto %%bash # Copy expression data table from GS to cluster gsutil cp gs://donnawerlingucsf/devEQTL/misc_files/counts.gene.cleaned-to-Donna-02092018-transpose.tsv counts.gene.cleaned-to-Donna-02092018-transpose.tsv ``` ```auto # Import sample annotation table with full expression data as pandas data frame exprdf = pd.read_csv('counts.gene.cleaned-to-Donna-02092018-transpose.tsv', sep='\t', index_col=0) ``` ```auto # Samples from VDS sampsInVds = vds.sample_ids ``` ```auto # Subset the expression data frame exprdf = exprdf.loc[sampsInVds,:] ``` ```auto # Calculate per-gene expression sums across this subsetted sample set exprsums = exprdf.sum(axis=0) ``` ```auto # Find non-expressed genes noExpr = exprsums[exprsums.le(0)] ``` ```auto # Remove these genes from the expression data exprdf = exprdf.drop(noExpr.index, axis=1) ``` ```auto # Write this df out to file exprdf.to_csv(path_or_buf='expressionFilt_176s_40707g-180411.txt', sep='\t') ``` ```auto %%bash # Transfer file to GS gsutil cp expressionFilt_176s_40707g-180411.txt gs://donnawerlingucsf/devEQTL/misc_files/ ``` ```auto %%capture # Re-load this expression data object as a keyTable exprdat = hc.import_table('gs://donnawerlingucsf/devEQTL/misc_files/expressionFilt_176s_40707g-180411.txt', impute=True).key_by('Sample') ``` ```auto # Add expression to sample annotations for this VDS vds = vds.annotate_samples_table(exprdat, root='sa.exprdat') ``` ```auto # Import sample annotation table with sample metadata information metadat = hc.import_table('gs://donnawerlingucsf/devEQTL/misc_files/eQTL_metadata_forDex-format-180316.txt', impute=True).key_by('Braincode') ``` ```auto # Add metadata to sample annotations for this VDS vds = vds.annotate_samples_table(metadat, root='sa.info') ``` ```auto # Add integer-coded variable for Sex vds = vds.annotate_samples_expr('sa.info.sexCode = if (sa.info.Sex == "f") 1 else if (sa.info.Sex == "m") 0 else NA: Int') ``` ```auto # Import gene locations file geneloc = hc.import_table('gs://donnawerlingucsf/devEQTL/misc_files/geneLocations.geneIdForHail.gencode.v21.txt', impute=True).key_by('geneId') ``` ```auto # Convert geneloc keyTable to pandas dataframe genedf = geneloc.to_pandas() ``` ```auto %%capture # Run the full loop, on all 40707 genes res_table = pd.DataFrame() for i in range (0, len(expr_genes)-1): # gene = '`' + expr_genes[i] + '`' testPheno = 'sa.exprdat.' + gene # chrVar = list(genedf.loc[genedf['geneId'] == gene]['chr'])[0] # startVar = list(genedf.loc[genedf['geneId'] == gene]['start'])[0] # endVar = list(genedf.loc[genedf['geneId'] == gene]['end'])[0] # Run linear regression for this gene vds_res = vds.filter_intervals(Interval(Locus(chrVar, startVar-1000000), Locus(chrVar, endVar+1000000)))\ .linreg(testPheno, covariates=['sa.info.Period','sa.info.sexCode']) # Pull out the results for this cis-eQTL analysis, filter to p<0.05, rename columns, convert to pandas data frame this_table = vds_res.variants_table().flatten()\ .select(['v','va.rsid','va.linreg.beta','va.linreg.se','va.linreg.tstat','va.linreg.pval'])\ .filter("`va.linreg.pval` > 0.05", keep=False)\ .rename(['Variant','rsID','Beta','SE','Tstat','P'])\ .annotate('Variant = str(Variant)') \ .to_pandas(expand=False) # Check - are there any significant results for this gene? nSig = this_table.size if nSig >= 1: # Add gene name to this table this_table['GeneId'] = gene.replace('`','') this_table[['EnsemblID','GeneSymbol']] = this_table['GeneId'].str.split('|',expand=True) # Split up the information in the variant column, re-order columns this_table[['Chr','Start','Ref','Alt']] = this_table['Variant'].str.split(':',expand=True) this_table = this_table[['GeneId','EnsemblID','GeneSymbol','Variant','Chr','Start','Ref','Alt','rsID','Beta','SE','Tstat','P']].head() # Add results to existing master results table res_table = pd.concat([res_table,this_table]) ``` ---
**Author:** ![dwerling](https://yyz2.discourse-cdn.com/flex036/user_avatar/discuss.hail.is/dwerling/32/138_2.png) [@dwerling](https://discuss.hail.is/u/dwerling)\ **Post date:** [April 11, 2018, 11:43pm UTC](https://discuss.hail.is/t/fast-linear-regression-for-multiple-phenotypes-eqtls/190/18 "2018-04-11T23:43:33Z")
And here is the error message I see, after the loop has processed about 9k of the 40k genes in the analysis: ```auto timeoutTraceback (most recent call last) in () 28 29 # Pull out the results for this cis-eQTL analysis, filter to p<0.05, rename columns, convert to pandas data frame ---> 30 this_table = vds_res.variants_table().flatten() .select(['v','va.rsid','va.linreg.beta','va.linreg.se','va.linreg.tstat','va.linreg.pval']) .filter("`va.linreg.pval` > 0.05", keep=False) .rename(['Variant','rsID','Beta','SE','Tstat','P']) .annotate('Variant = str(Variant)') .to_pandas(expand=False) 31 32 # Check - are there any significant results for this gene? in to_pandas(self, expand, flatten) /home/hail/hail-0.1-74bf1ebb3edc.zip/hail/java.pyc in handle_py4j(func, *args, **kwargs) 113 def handle_py4j(func, *args, **kwargs): 114 try: --> 115 r = func(*args, **kwargs) 116 except py4j.protocol.Py4JJavaError as e: 117 tpl = Env.jutils().handleForPython(e.java_exception) in to_pandas(self, expand, flatten) /home/hail/hail-0.1-74bf1ebb3edc.zip/hail/typecheck/check.pyc in _typecheck(f, *args, **kwargs) 243 def _typecheck(f, *args, **kwargs): 244 check_all(f, args, kwargs, checkers, is_method=True) --> 245 return f(*args, **kwargs) 246 247 return decorator(_typecheck) /home/hail/hail-0.1-74bf1ebb3edc.zip/hail/keytable.pyc in to_pandas(self, expand, flatten) 642 """ 643 --> 644 return self.to_dataframe(expand, flatten).toPandas() 645 646 @handle_py4j /usr/lib/spark/python/pyspark/sql/dataframe.pyc in toPandas(self) 1440 """ 1441 import pandas as pd -> 1442 return pd.DataFrame.from_records(self.collect(), columns=self.columns) 1443 1444 ########################################################################################## /usr/lib/spark/python/pyspark/sql/dataframe.pyc in collect(self) 309 with SCCallSiteSync(self._sc) as css: 310 port = self._jdf.collectToPython() --> 311 return list(_load_from_socket(port, BatchedSerializer(PickleSerializer()))) 312 313 @ignore_unicode_prefix /usr/lib/spark/python/pyspark/rdd.pyc in _load_from_socket(port, serializer) 140 try: 141 rf = sock.makefile("rb", 65536) --> 142 for item in serializer.load_stream(rf): 143 yield item 144 finally: /usr/lib/spark/python/pyspark/serializers.pyc in load_stream(self, stream) 137 while True: 138 try: --> 139 yield self._read_with_length(stream) 140 except EOFError: 141 return /usr/lib/spark/python/pyspark/serializers.pyc in _read_with_length(self, stream) 154 155 def _read_with_length(self, stream): --> 156 length = read_int(stream) 157 if length == SpecialLengths.END_OF_DATA_SECTION: 158 raise EOFError /usr/lib/spark/python/pyspark/serializers.pyc in read_int(stream) 541 542 def read_int(stream): --> 543 length = stream.read(4) 544 if not length: 545 raise EOFError /usr/lib/python2.7/socket.pyc in read(self, size) 378 # fragmentation issues on many platforms. 379 try: --> 380 data = self._sock.recv(left) 381 except error, e: 382 if e.args[0] == EINTR: timeout: timed out ``` I would also note that this loop actually runs to completion on my laptop, in about 3 hours (this is on a variant-downsampled VDS). The timeout issue only happens when I try to run the same analysis on a dataproc cluster. ---
**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:** [April 12, 2018, 12:25am UTC](https://discuss.hail.is/t/fast-linear-regression-for-multiple-phenotypes-eqtls/190/19 "2018-04-12T00:25:59Z")
Thanks for the detail! That’ll make this a lot easier to figure out. This looks similar: [https://github.com/databricks/spark-sklearn/issues/49](https://github.com/databricks/spark-sklearn/issues/49) I think it might also make things a bit easier if you export the linreg results to a file instead of bringing back to pandas, and then concatenate everything later. That way, if things fail again, you can restart from the point where it failed! I’ll try to find time for a reply commensurate with your detailed report, but it might not be until Friday. ---
**Author:** ![nnekaede](https://yyz2.discourse-cdn.com/flex036/user_avatar/discuss.hail.is/nnekaede/32/155_2.png) [@nnekaede](https://discuss.hail.is/u/nnekaede)\ **Post date:** [July 13, 2018, 7:24pm UTC](https://discuss.hail.is/t/fast-linear-regression-for-multiple-phenotypes-eqtls/190/20 "2018-07-13T19:24:34Z")
Is this function available in Hail 0.2? [Next page](https://discuss.hail.is/t/fast-linear-regression-for-multiple-phenotypes-eqtls/190.md?page=2)