Error message: system is computationally singular

[evatosco]

Hi!

I'm trying to run BHR in univariate mode and grouping by genes, from a VCF file containing my filtered interest variants (around 3k). After adjusting the code and trying to run it, I encountered this error message:

Burden Heritability Regression
Daniel Weiner, Ajay Nadig, and Luke O'Connor, 2023
Running BHR at 2023-09-04 15:15:21.595497
For trait 1, MAF ranges from 0.00061 to 0.49
MAF ranges over more than two orders of magnitude for trait 1. We recommend splitting into finer MAF bins. See documentation
Lambda GC:  1.524
...seems reasonable
39 genes in the BHR regression, 1 of which are significant fixed effects
...seems reasonable

Error in solve.default(matrix(rowSums(jackknife_denominators), ncol = sqrt(length(rowSums(jackknife_denominators))),  : 
  system is computationally singular: reciprocal condition number = 4.37745e-18

I've read here that it could be a limitation in computation methods, but I have no idea about the detail of this or how to solve it.

Would anyone have an insight about this? I cannot share the input data for privacy reasons but would be happy to answer any questions or suggestions.

My txt input file (from the VCF file) to the BHR() function has this structure:

gene    AF      beta    gene_position chromosome N  phenotype_key
ENSG001 0.0006082725 -0.703667280      45810531          1 822          pass
ENSG001 0.0006082725 -0.703667280      45810560          1 822          pass
ENSG001 0.0006082725  1.422856632      45810646          1 822          pass
ENSG001 0.0006082725  1.422856632      45810659          1 822          pass
ENSG001 0.0006082725  1.422856632      45810774          1 822          pass
ENSG001 0.0006082725 -0.703667280      45810797          1 822          pass
ENSG001 0.0006082725  1.422856632      45810808          1 822          pass

Thanks! Eva

[danjweiner]

Hi Eva,

Thank you for your question, and apologies for the delayed response.

The computationally singular error is appearing because there are too few genes in the BHR regression you are implementing. Specifically, there are 39 genes in your regression, which is fewer than the default number of jackknife blocks (n = 100) in the model; this will lead to some blocks having 0 genes in them, which will raise this error.

A few questions and comments:

• In your input summary statistics, are there 39 genes or are some genes being lost between your input statistics and the model run? BHR is a regression across genes, and it will perform poorly when the number of genes is small.
• You have included in a single BHR run variants across a wide range of allele frequencies; linkage disequilibrium between variants will inflate burden heritability estimates. We recommend stratifying variants by frequency and running regressions separately across frequency groups (see manuscript for more details). Specifically, analysis of rare variants (in our paper, AF < 1e-3) can help reduce the influence of LD on burden heritability estimates.
• It also appears you are running BHR with summary statistics from a small n association study (822, not sure if this a quantitative trait or case/control). BHR will likely be underpowered; in our analysis of neuropsychiatric traits in the paper, we ran into power considerations around 10k cases.

Please let us know if you have additional questions!