Genome-wide meta-analysis of insomnia in over 2.3 million individuals implicates involvement of specific biological pathways through gene-prioritization

SNP and gene-based findings from the meta-analysis

The meta-analysis yielded 51,876 genome-wide significant SNPs, residing in 554 distinct loci containing 791 independent lead SNPs (r²<0.1; Methods, Extended Data Fig. 2-3, Supplementary Table 3-4 and Supplementary Note 3). Of the top 538 SNPs (i.e. SNPs with the lowest P-value in each locus), 97.1% showed concordance in direction of effect between the two cohorts. Out of 554 loci, 11 loci were genome-wide significant in the UKB, and 419 also in the 23andMe cohort-specific GWASs, while 9 loci were identified in both cohorts (Supplementary Table 4). The total summed length of risk loci was 145.2Mb, which is 4.9% of the genomic regions containing known SNPs in the entire genome. Of the 554 loci, 190 loci overlapped with previously identified risk loci^5–7,11 and 364 were novel (Methods, Fig. 1a and Supplementary Table 4). Of the loci previously reported in Hammerschlag et al.¹¹, Lane et al. (2017)⁵, Lane et al. (2019)⁷, and Jansen et al.⁶, 1/2, 3/5, 25/57 and 18/202 loci were no longer significant in the current meta-analyses, respectively (Supplementary Table 4, Supplementary Note 4 and Extended Data Fig. 4). We show that the number of risk loci increases almost linearly as a function of sample size (Fig. 1b), and that both newly identified risk loci and un-replicated risk loci from the previous GWASs showed significantly higher P-values compared to other risk loci, as expected (Supplementary Note 5 and Extended Data Fig. 5).

• Download figure
• Open in new tab

Figure 2. Tissues, brain regions, cell types and gene-sets associated with insomnia based on MAGMA analyses.

P-values of significantly associated tissues, brain regions, cell types and gene-sets. Confounders after pairwise conditional analyses per dataset are not displayed (full results are available in Supplementary Tables 13-21). Horizontal solid lines above bars represent clusters of items which are either collinear or jointly associated with insomnia. Bars are colored by datasets. A letter in the parenthesis for brain regions from the AHBA high dataset refers to parent brain regions from the AHBA low dataset (F: frontal lobe, O: occipital lobe, P: parietal lobe and T: temporal lobe). Labels in parenthesis for cell types from DropViz dataset refer to the cluster label from the original study. BP: biological process, CC: cellular components. The dashed line indicates the Bonferroni corrected threshold for statistical significance (p=0.05/5974).

• Download figure
• Open in new tab

Figure 3. Schematic overview of gene prioritization strategy from risk loci.

(a) a conceptual illustration of gene-prioritization based on multi-loci information. Rectangles represent genes colored by prioritization status, HC-1: high confidence, single gene from a locus; HC-m: high confidence in locus with multiple HC-m genes. Genes colored grey are located within or close to insomnia risk loci but were not implicated either due to a lack of functional evidence or because other genes in the same locus had higher priority. Genes that do not have a direct protein-protein interaction to HC-1 genes are filtered out (indicated by a cross). (b) a flowchart of gene prioritization using credible SNPs and genome-wide significant (GWS) SNPs. Solid boxes represent inputs or outputs, dashed grey boxes represent conditions, green arrows represent paths where an originated condition is positive, orange arrows represent paths when an originated condition is negative, grey arrows represent paths which do not have any other option.

• Download figure
• Open in new tab

Figure 4. Tissues and gene-sets associated with the high-confidence prioritized genes.

P-values of significantly associated tissues, brain regions, cell types and gene-sets. Confounders after pairwise conditional analyses per dataset are not displayed. Horizontal solid lines above bars represent cluster of items which are either collinear or jointly associated with insomnia. Stars represent items that also showed significant association with insomnia by using the full GWAS results. Bars are colored by datasets. BP: biological process, CC: cellular components, MF: molecular function.

The 51,876 genome-wide significant SNPs were most abundant in intronic regions (47.7%) and significantly enriched compared to all analyzed SNPs (enrichment (E)=1.28; p<1e-323, two-sided Fisher’s exact test; Methods). Intergenic regions (34.4%) were the second most common, although these showed a significant depletion compared to all analyzed SNPs (E=0.75, p<1e-323; Fig. 1c and Supplementary Table 5). GWS SNPs in exonic regions (0.73%) showed a significant depletion (E=0.86, p=3.6e-3) compared to all analyzed SNPs, and included 168 nonsynonymous, 1 stop gain, 192 synonymous and 1 non-frameshift substitution. SNPs in 3’ untranslated regions (UTR; 1.1%), representing significant enrichment (E=1.30, p=6.9e-10), while proportions of SNPs in other flanking regions did not differ from all analyzed SNPs (Fig. 1c and Supplementary Table 5). Stratified heritability analyses¹⁴ with 28 annotations (Methods), showed that SNP heritability was most strongly enriched in conserved regions, followed by multiple chromatin modification markers, consistent with the previous meta-analysis⁶ (Fig. 1d and Supplementary Table 6). Next, we performed a gene-based test using MAGMA gene analysis by assigning SNPs to genes with 2 kb upstream and 1 kb downstream windows. Of 19,751 protein-coding genes analyzed, 1,429 genes reached genome-wide significance (0.05/19751=2.53e-6; Fig. 1a and Supplementary Table 7). The most significant gene was PTPRD (p=7.2e-37) which has been associated previously with insomnia ⁶, restless leg syndrome¹⁵, type 2 diabetes¹⁶, and coronary artery disease¹⁷. The second-most significant gene was LSAMP (p=2.8e-36) which was not significant in the previous insomnia GWAS but has been associated with major depressive disorder (MDD)¹⁸ and suicidal behavior¹⁹, which are highly genetically correlated to insomnia. The most significantly associated genes from the previous insomnia GWASs, MEIS1¹¹ and BTBD9⁶, were also supported in the current study (p=1.2e-14 and 4.8e-24, respectively). Risk loci and gene analyses for sex-specific meta-analyses are summarized in Supplementary Note 6, Supplementary Table 7-11 and Extended Data Fig. 6.

The 51,876 GWS SNPs were mapped to 3,526 genes (of which 1,455 are located within the risk loci) using positional, eQTL and chromatin interaction mapping strategies²⁰ (Methods and Supplementary Table 12). Together with genes significantly associated with insomnia based on formal gene-based analysis in MAGMA, 3,898 unique genes were implicated. We observed significant genetic correlations with 350 traits out of 551 traits tested, including multiple cardiovascular, metabolic, and psychiatric traits, in agreement with previous reports⁶ (Methods, Supplementary Note 7, Supplementary Table 13 and Extended Data Fig. 7). Of 554 insomnia risk loci, 282 loci were colocalized with one of the 350 traits indicating shared causal variants between the traits (Supplementary Note 7). In addition, a clustering of these 282 loci based on the colocalization pattern across 350 traits suggested presence of locus heterogeneity where we observed distinct clusters of loci: one was mainly colocalized with metabolic traits and the other was mainly colocalized with psychiatric traits (Supplementary Note 7, Supplementary Table 14-18 and Extended Data Fig. 8). We performed extensive post-GWAS analyses to test for convergence of genetic association signals in tissue types, brain regions, cell types and biological pathways associated with insomnia, based on the total genome-wide distribution of genetic associations, weighted by statistical significance using MAGMA gene-set analyses²¹ (Methods and Supplementary Note 8). Based on the full GWAS distribution, we found evidence for enrichment in brain tissue (specifically the cerebral cortex), neurons in four brain specific areas (lateral geniculate nucleus, habenula, ventral pallidum, and anterior pretectal nucleus), GABAergic neuron and biological pathways involved in synaptic functions, behavior and neuron differentiation (Fig. 2, Supplementary Note 8, Extended Data Fig. 9 and Supplementary Tables 19-30).

Prioritization of high confidence genes from 554 risk loci using multi-locus information

Due to LD, many non-causal SNPs will show a statistical association with a trait simply because they are correlated with causal SNPs due to LD. We call these non-causal significant SNPs ‘LD by-products’. Since we do not know which SNPs are the causals ones and which are the LD by-products, post-GWAS annotation will provide information on many SNPs that are likely not causally related to the trait, and therefore conventional post-GWAS analysis testing for convergence may still contain considerable noise from these LD by-products. Multiple studies have already proposed to conduct fine-mapping per locus and prioritize credible SNPs and genes from each locus before testing for convergence^22–25. Here we propose to prioritize genes additionally based on cross-locus connections.

We assume the following: (i) Credible SNPs can be indicated using in silico fine-mapping strategies; (ii) SNPs that have a structural or regulatory effect on a gene product are more likely to be causal than SNPs that have no such (known) effect; (iii) genes that are implicated as the only gene in a locus are likely to be the gene responsible for the statistical result at the locus, and are therefore likely a true causal gene, and (iv) if insomnia is influenced by hundreds of genes, and at least part of those are functionally related (Fig. 3).

We first defined credible SNPs by performing statistical fine-mapping for each of the 554 insomnia risk loci, using FINEMAP²³ (Methods). For each locus, the 95% credible sets were extracted, resulting in a total of 26,016 unique SNPs (see Supplementary Note 9 and Extended Data Fig. 10 for detailed results). More than 94.5% of the SNPs in 95% credible sets had a posterior inclusion probability (PIP) ≤0.1, suggesting that those SNPs were ‘unsolved’ by FINEMAP, which can be due to high LD and small effect sizes. We retained only the 1,423 credible SNPs with PIP>0.1 distributed over 429 loci. Association statistics and functional annotation of these credible SNPs are provided in Supplementary Table 31. Credible SNPs were then mapped to genes if they were deleterious coding SNPs (non-synonymous, stop-gain, stop-loss or splicing SNPs) (Supplementary Table 32) or if they were eQTLs from GTEx v8²⁶, PsychENCODE²⁷ and eQTLGen²⁸ resources that were colocalized with insomnia GWAS summary statistics (Methods and Supplementary Table 33). This resulted in labeling 314 genes from 178 loci as having a credible SNP that was either nonsynonymous coding or co-localized with eQTLs (Supplementary Table 34). For 376 of the remaining insomnia loci, no genes were present with the above criteria. This was either due to a lack of credible SNPs or a lack of credible SNPs that were deleterious coding or were colocalized eQTLs. For these loci, we used GWS SNPs instead of credible SNPs and again evaluated whether they were deleterious coding SNPs or colocalized with eQTLs. This resulted in an additional 257 genes from 103 loci (Supplementary Table 34). Together, these 571 genes from 281 loci were called high confidence (HC) genes. A single locus could contain multiple HC genes, but some were mapped to a single gene. In the case of a single HC gene, that gene was considered the only likely causal candidate from the locus. We labeled these single HC genes as HC-1 (with the highest confidence) and the remaining genes as HC-m (i.e., HC genes but with multiple genes in the same locus and therefore unclear whether they are all causal, whether only one is causal, or whether a subset is causal; Fig. 3 and Methods). We obtained 164 HC-1 genes from 166 loci (2 genes were mapped from two loci due to distal eQTLs) and 407 HC-m genes from 116 loci (Supplementary Table 34).

Assuming that the HC-1 genes are the most likely causal genes, and that functional relationships exist among the set of causal genes, we then used the set of HC-1 genes to select the most likely causal genes from the set of HC-m genes. To assess functional relations, we used protein-protein interactions (PPI) using the InWeb database²⁹, and identified HC-m genes that have a direct interaction with HC-1 genes (Methods). We choose to use PPI for this, and not for example co-expression, because we aimed to identify genes whose products are known to form a protein complex, increasing the likelihood that they have a common function. Alternate measures of connections between genes (such as spatial or temporal co-expression) can of course also be used in this step. Of 407 HC-m genes, 125 genes from 74 loci were found to have a direct PPI with an HC-1 gene. We then defined 289 genes (164 HC-1 and 125 HC-m genes connected to HC-1) from 239 loci as the ‘high confidence prioritized’ (HCP) genes (Supplementary Table 34). The HCP genes are the genes we believe to be most likely functionally associated to insomnia based on the associated loci. Out of 239 loci, 202 were linked to single HCP genes and the maximum number of genes from a single locus was 5. We consider these 239 loci to be ‘resolved’ and the remaining 315 loci to be unsolved due to lack of biological evidence or information at this moment. Although we still cannot identify genes from 56.7% of the insomnia loci, we assume that these loci are randomly distributed, and thus we expect the HCP genes to explain some of the underlying biological mechanisms of insomnia as the HCP genes are more refined and likely to contain fewer false positives compared to 3,526 genes mapped by all GWS SNPs.

The HCP genes include NEGR1 which has been reported to be associated with multiple traits such as body mass index³⁰, MDD³¹, cognitive traits^32,33 as well as in the previous insomnia GWAS⁶ and its effect on neuronal growth and behavior has been suggested^34,35. The HCP genes also include EP300 known to be involved in circadian rhythm³⁶ which regulates sleep timing (see Supplementary Table 34 for the full list of the HCP genes).

To look for convergence in biological functions of the HCP genes, we assessed enrichment in tissues, brain regions, cell types and biological pathways (Methods).

The HCP genes showed significant joint associations with frontal cortex BA9 and anterior cingulate cortex BA24 (Fig. 3 and Supplementary Table 35-36). These tissue enrichments seen with HCP genes were also detected when using the full GWAS results (Supplementary Table 37-39 and Supplementary Note 8).

On the other hand, some tissue enrichments that were significant when using the full GWAS, are no longer detected when using the HCP gene list. This is because genes driving these associations in the full GWAS were not prioritized, which suggests that these initial enrichments were less likely to be based on true causal genes. (Supplementary Note 10 and Extended Data Fig. 11).

From the gene-set enrichment analyses, the HCP genes showed significant enrichment in 18 gene-sets clustered into four independent groups (Fig. 4, Extended Data Fig. 12 and Supplementary Table 40). The most significantly enriched gene-sets from each cluster were ‘modulation of chemical synaptic transmission’ (SynGO:BP), ‘neuron differentiation’ (GO:BP), ‘regulation of trans synaptic signaling’ (GO:BP) and ‘inclusion body’ (GO:CC). We identified not completely overlapping and more gene-sets when using the prioritized genes than when using the full GWAS results. This may indicate that, by filtering LD by-products genes to reduce noise in associated genes, we were able to identify gene-sets specific to the most likely causal genes prioritized from insomnia associated loci, which are missed by using full GWAS results.

In this study, we performed a GWAS meta-analysis for insomnia including over 2.3 million subjects. We identified 554 insomnia risk loci, almost doubling the number identified by the previous largest study with 1.3 million subjects⁶. We demonstrated a novel strategy using known biological functions of SNPs and multi-locus functional relations of genes to prioritize the most likely causal genes, and based post-GWAS analyses for convergence on these genes. Applying this strategy, we identified 289 HCP genes from 239 loci and compared associated tissue and cell types, as well as gene-sets based both on the set of prioritized genes and all genes implicated in the GWAS. We found that the former is less likely to contain LD by-products, as it provided more specific results, pointing towards specific synaptic functions in specific brain regions. Of course, the ultimate test of whether the HCP genes are actually causally involved still lies in functional follow-up experiments.

There are several limitations of the current strategy. First, the proposed strategy using multi-locus information is only feasible for polygenic traits with a reasonable number of independent risk loci identified by GWAS and a reasonable number of loci with single implicated genes. Second, the prioritization procedure depends on the availability and accuracy of functional annotations of SNPs and genes. For example, we defined an HC-1 gene as the only gene from a single locus with high confidence biological evidence.

However, in the future, more SNPs may be known to be deleterious (by increasing accessible (rare) SNPs in GWAS) or colocalized eQTLs (by increasing statistical power to detect eQTLs and availability of them in specific cell types), which may change the current results and may allow us to identify additional high confidence genes from loci. Third, cross-locus linking of genes depends on the availability and reliability of biological information (PPI, co-expression networks, or any other gene-correlation matrix deemed relevant), which is currently not abundantly available and still imperfect. We do believe using cross-locus information greatly aids in making sense of the multitude of associated genes, and the current study shows that this strategy implies a role for more specific biological functions in insomnia. In conclusion, the current study provides further insights into the possible genetic pathways contributing to insomnia using a novel strategy in prioritizing genes and generating hypotheses that can be tested in functional experiments.

Sample cohorts

Phenotype assessment

Genome-wide association analysis

GWAS meta-analysis

A meta-analysis on the GWAS summary statistics from UKB and 23andMe cohorts was performed using METAL software¹², based on the fixed-effect model with SNP P-values weighted by sample size. The meta-analysis was performed for the sex-combined and each sex-specific GWASs separately. We used a unique marker ID consisting of chromosome, position, and alphabetically ordered alleles to match SNPs between cohorts. Since indels were coded as I/D in the 23andMe cohort, exact alleles were assigned for indels based on information from the UKB cohort. We assigned alleles to indels in the 23andMe only when there were bi-allelic indels on the same position in the UKB cohort. The number of analyzed SNPs and sample sizes are summarized in Supplementary Table 1.

We converted Z-statistics to standardized effect sizes (and their standard error) as a function of MAF and sample size as below⁴⁹. where p is MAF, n is sample size for a given SNP.

SNP heritability and stratified heritability estimate

SNP heritability was estimated for sex-combined and sex-specific GWASs on UKB, 23andMe, and meta-analyses, using LD score regression (LDSC)¹³. Pre-computed LD scores for 1000 Genome Phase 1 European subjects were downloaded from https://data.broadinstitute.org/alkesgroup/LDSCORE/. The analyses were limited to

HapMap3 SNPs and the MHC region (chr6:26Mb-34Mb) was excluded. Additionally, SNPs with chi-square statistics >80 were excluded. To compute SNP heritability on the liability scale, we provided a population prevalence of 30%¹. Since the sample size of this study is large, λ_GC was scaled for 1000 cases and 1000 controls as λ₁₀₀₀ = 1 + (λ_GC - 1) x (1/Ncases + 1/Ncontrols) x 500.

To test whether SNP heritability was enriched in a specific category of functional annotations, we partitioned SNP heritability for binary SNP annotations¹⁴. The enrichment was computed as the proportion of SNP heritability explained by the SNPs with the annotations divided by the proportion of SNPs with the annotations. We obtained 28 functional annotations from https://data.broadinstitute.org/alkesgroup/LDSCORE/.

Polygenic risk scoring

To evaluate phenotypic variance explained by our insomnia meta-analysis, polygenic score (PGS) was computed based on the SNP effect sizes using PRSice v2.2.1⁵⁰ (http://www.prsice.info/) with default parameters. PRSice performs clumping of SNPS at r²=0.1. We set P-value thresholds of input SNPs at p<1.0, 0.5, 0.05, 0.01, 5-e3, 1e-3 and 1e-5. We computed PGS for three randomly selected sets of 10,000 UKB subjects, and summary statistics were re-calculated excluding those 10,000 subjects each time for the UKB cohort. We then meta-analyzed with the 23andMe GWAS which we used as training data. To compare the difference in the size of the target subjects, we also computed PGS for 3,000 subjects, 3 times by the same procedure. Chromosome X was excluded from the analyses. To evaluate the predictive power of the insomnia GWAS meta-analysis, we computed PGS for an independent cohort, the Million Veteran Program (MVP). The MVP is a research program funded by the US Department of Veterans Affairs to learn how genes, lifestyle, and military exposures affect health and illness⁵¹. In the MVP cohort, insomnia cases were defined as individuals answering “Yes” to any of the following three questions : “Trouble falling asleep when you first go to bed”, “Waking up during the night and not easily going back to sleep, or “Waking up in the morning earlier than planned or desired”. To be defined as a case, individuals also had to respond “Yes” to the following: “Feeling unsatisfied or not rested by your night’s sleep”. Controls were defined as not responding “Yes” to any of the first three questions and responding “No’ to the final question. This resulted in 45,355 cases and 138,589 controls. We then computed PGS using PRSice in the same way as described above, except the summary statistics of the full meta-analysis were used as training set this time.

Definition of risk loci

Genomic risk loci were defined within FUMA (https://fuma.ctglab.nl) as described previously²⁰, by first clumping SNPs with p<5e-8 at r²=0.6 using SNPs with p<1e-5 to define independent significant SNPs. Those independent significant SNPs were further clumped at r²=0.1 to define lead SNPs. The independent significant SNPs that were in LD with the same lead SNPs (r²>0.1) and LD blocks closer than 250kb were merged into a single locus. Each locus is represented by the most significant SNP (top SNP). A risk locus can contain multiple lead and independent significant SNPs. We manually excluded suspicious loci (e.g., a single SNP reaching genome-wide significance without any SNPs with p<1e-5 with r²>0.6) by examining LocusZoom plots (see Supplementary Note 2 for details).

Comparison with previously identified insomnia risk loci

Previously identified insomnia risk loci were obtained from 4 studies; 2 loci from Hammerschlag et al.¹¹, 5 loci from Lane et al. (2017)⁵, 57 loci from Lane et al. (2019)⁷ and 202 loci from Jansen et al.⁶. We then overlapped previously identified loci with 554 loci identified from the sex-combined meta-analysis in this study. We considered risk loci which did not overlap with any of previously identified loci as novel findings.

Comparison of insomnia GWASs with different sample sizes

Insomnia GWAS summary statistics for 6 different sample sizes were obtained as below. UKB1 (N=113,006): The summary statistics were obtained from the study of Hammerschlag et al.¹¹ (https://ctg.cncr.nl/software/summary_statistics).

UKB2 (N=386,988): UKB only summary statistics from the current study. 23andMe1 (N=944,477): Only 23andMe cohort from the study of Jansen et al.⁶ (summary statistics is not publicly available, only the number of risk loci was obtained).

Meta1 (N=1,331,010): The summary statistics of meta-analysis from the study of Jansen et al.⁶ (summary statistics is not publicly available, only lead SNPs were obtained).

23andMe2 (N=1,978,022): 23andMe only summary statistics from the current study (not publicly available).

Meta2 (N=2,365,010): The summary statistics of the main meta-analysis from the current study (not publicly available).

For each GWAS, risk loci were defined as described above. For 23andMe1 and Meta1, the number of risk loci were obtained from the previous study⁶ as the summary statistics are not available. For comparison of P-values of the same locus across GWASs with different sizes, we obtained lead SNPs for Meta1 from the previous study⁶ while 23andMe1 was excluded from the analysis.

To obtain matched SNPs across GWAS to compare effect sizes, we first extracted genome-wide significant SNPs for all risk loci from 5 GWASs (23andMe1 was excluded from the analysis). For Meta1, the analysis was limited to lead SNPs provided by the previous study⁶, therefore not all GWS SNPs in the risk loci were assessed. We then selected SNPs that reached genome-wide significance in at least 3 GWASs out of 5. From each locus, the single SNP with the minimum P-value was further selected (when there are multiple) for the analyses. Standardized effect sizes of those SNPs were computed based on summary statistics from each GWAS as described above. When Z-score were not available, P-values were converted to Z-score (two-sided).

MAGMA gene, gene-property and gene-set analyses

Gene-base testing was performed using MAGMA v1.07²¹ to obtain gene P-values, using summary statistics of sex-combined meta-analysis. From 20,260 protein-coding genes, SNPs were assigned to one of the 19,751 genes within 2kb up- and 1kb downstream windows, based on the location obtained from Ensembl v92 GRCh37 using BioMart. We used the SNP-wise mean model and randomly selected 10,000 unrelated European subjects from the UKB cohort were used as a reference panel.

For tissue specificity analyses, we obtained RNA-seq data for 54 tissue types from GTEx v8²⁶. Read per kilobase per million (RPKM) was log-2 transformed with pseudo count 1 followed by winsorization at 50, and average per tissue type was computed for each gene. In a gene-property analysis, the average across 54 tissue types was conditioned and the one-sided test was performed to identify positive associations of tissue specific gene expression with insomnia.

For gene expression of specific brain regions, we obtained normalized microarray data for 3,702 samples from 6 healthy donors⁵², from Allen Human Brain Atlas (AHBA; http://human.brain-map.org/static/download). From 58,692 probes, 31,098 with missing values in less than 20% of the samples were first extracted. When there were multiple probes per gene, genes with the highest variance across 3,702 samples were selected, resulting in 17,916 unique genes. Of those, 13,943 genes were mapped to unique Ensembl gene ID (v92

GRCh37). Every 3,702 samples were assigned to the structural ID of the brain where there were multiple layers of hierarchical structure⁵². We assigned the annotation of brain regions at 4th and 5th layers of the hierarchical structure where the top of the tree (‘brain’) was considered as layer 0. We limited consideration to brain regions with at least 5 samples, resulting in 54 and 106 brain regions for layer 4 and 5, respectively. For each brain region, the average expression values were computed for each gene. We then performed gene-property analysis conditioning on average across brain regions and a one-sided test was performed to identify positive associations of brain region specific gene expression with insomnia.

For cell type specificity analysis, we used 5 datasets from Linnarsson’s group that were used in the previous insomnia meta-analysis study: mouse samples from cortex and hippocampus (GSE60361, level 2 neurons)⁵³, hypothalamus (GSE74672, level 2 neurons)⁵⁴, oligodendrocytes (GSE75330)⁵⁵, midbrain (GSE76381)⁵⁶ and striatum (GSE97478)⁵⁷. We additionally tested one of the most comprehensive scRNA-seq datasets for mouse brain, DropViz⁵⁸ (http://dropviz.org/, 565 sub-clusters from 9 brain regions). In total, we obtain 728 cell types from 6 datasets. Each dataset was pre-processed as described in the previous study⁵⁹. We performed the 3-step workflow (per dataset analysis, within-dataset and cross-dataset conditional analyses) as described previously⁵⁹ to determine significant associations supported by multiple independent datasets..

Gene-sets for Gene Ontology (GO) terms and canonical pathways were obtained from MsigDB v6.2⁶⁰ (http://software.broadinstitute.org/gsea/msigdb/index.jsp). We performed gene-set analyses for 5,033 gene-sets with at least 20 genes.

We evaluated all 5,974 tested items against the Bonferroni corrected threshold (0.05/5974=8.4e-6).

Fine-mapping and credible SNPs

We performed fine-mapping of 554 risk loci using FINEMAP with the shotgun stochastic search algorithm²³ (http://www.christianbenner.com/). For each risk locus, SNPs within 50kb from the top SNP (with the minimum P-value) or locus boundary, whichever larger, and with p<0.05 were used for fine-mapping. The pairwise LD matrix of SNPs was estimated based on randomly selected 100k unrelated European individuals from UKB using LDstore⁶¹ (http://www.christianbenner.com/). The maximum number of causal SNPs (k) was set to 10.

An exception was made for locus #4 which contained 7 SNPs, where k was set to 5. FINEMAP outputs a set of models (all possible combination of k causal SNPs in a locus) with posterior probability (PP) of being a causal model. A 95% credible set was defined by taking models from the highest PP until the cumulative sum of PP reached 0.95 for each locus. The 95% credible set SNPs were, therefore, defined by taking unique SNPs from 95% credible sets. In addition, a posterior inclusion probability (PIP) was calculated for each SNP as the sum of PPs of credible sets that contain that SNP. In this study, we only used credible SNPs with PIP>0.1, however, results of all credible set SNPs are described in Supplementary Note 7 and Extended Data Fig. 8.

Annotation of SNPs with FUMA

All GWS SNPs were annotated with their functions using FUMA. Functional consequences of SNPs on genes were obtained by performing ANNOVAR⁶² gene-based annotation using Ensembl genes. Enrichments of GWS SNPs in each annotation were tested by Fisher’s exact test (two-sided) by compared in with the annotations of all SNPs analyzed in the meta-analysis. To determine the deleteriousness of the SNPs, CADD v1.4⁶³ score was annotated for each SNPs. In addition, RegulomeDB score (categorical score indicating how likely the SNPs are involved in regulatory elements)⁶⁴ and 15-core chromatin state for 127 tissue/cell types obtained from Roadmap^65,66 were annotated.

eQTL colocalization

We colocalized insomnia summary statistics from the sex-combined meta-analysis within 554 risk loci with eQTL summary statistics using the coloc package in R⁶⁷. We only tested genes whose significant eQTLs were overlapping with at least one of the GWS SNPs in the insomnia GWAS, and colocalization was performed for each gene for each of 51 eQTL datasets (i.e. 49 tissues from GTEx v8²⁶, meta-analysis of blood samples from eQTLGen consortium²⁸ and prefrontal cortex from PsychENCODE²⁷). For each colocalization, we extracted SNPs available in both eQTL summary statistics of a testing gene and within 10kb of insomnia risk loci. We then performed coloc.abf function. We did not perform colocalization when there were <10 SNPs overlapping between insomnia and eQTL summary statistics.

The coloc.abf function assumes a single causal SNP for each trait and estimates the posterior probability of the following 5 scenarios for each testing region; H₀: neither trait has a genetic association, H₁: only trait 1 has a genetic association, H₂: only trait 2 has a genetic association, H₃: both trait 1 and 2 are associated but with different causal SNPs and H₄: both trait 1 and 2 are associated with the same single causal SNP. In the case of our study, trait 1 is insomnia and trait 2 is expression of a tested gene. Since we limited the analyses to genes where there was at least one overlap of significant SNPs with insomnia GWS SNPs, this discards scenarios H₀ to H₂, and we are only interested in whether H₄ is most likely. We therefore defined eQTLs of a testing gene as colocalized with insomnia summary statistics when H₄ is >0.9. It is possible that genomic regions outside of the pre-defined risk loci can also be colocalized with eQTLs. However, we limited the analyses to within the risk loci in this study, since the primary aim was to prioritize genes linked from the insomnia risk loci.

Gene mapping with FUMA

We used FUMA to map SNPs to genes with 3 criteria: positional, eQTL and chromatin interaction mappings²⁰.

Tissue and cell type association, and gene-set enrichment tests with prioritized genes

Associations of tissue and cell type specific gene expression with prioritized genes were tested with a linear regression model using lm function in R. We defined the model as follows to correct for the average expression across tissues or cell types within a dataset and gene size: where G is a binary status reflecting whether the gene was prioritized (1) or not (0), E is a tissue/cell type specific expression value, A is the average expression of the gene across all available tissues/cell types in a dataset, S is a gene length. We performed a one-sided test (β_G>0) which evaluates how well the prioritized gene status predicts the specificity of gene expression in a testing tissue/cell type. For tissue/cell types significantly associated with the prioritized genes, we performed conditional analyses as below and performed a one-sided test (β_{G_1}>0 and β_{G_2}>0). For gene-set enrichment analyses, a one-sided hypergeometric test (greater) was performed. We tested the same datasets as tested in the MAGMA gene-set analysis (53 tissues, 54 and 106 brain regions, 728 cell types, and 5,033 gene-sets). As is done in MAGMA v1.07, gene-property (gene expression value) was truncated when the value is above or below 5 standard error from the mean. The analyses were limited to genes available in each dataset out of 20,260 protein-coding genes based on Ensembl v92 GRCh 37.

We performed Bonferroni correction across all 5,974 tested items (0.05/5974=8.4e-6).

Conditional analysis for associated tissues, cell types and gene-sets using the full GWAS results and the set of prioritized genes

To identify independent association signals in gene-set or gene-property associations, we performed conditional analyses for sets that were significantly associated with insomnia per dataset (i.e., GTEx, AHBA low/high resolution, scRNA and gene-sets) in both the analyses based on the full genome-wide distribution (using linear regression as implemented in MAGMA) and the set of prioritized genes (using linear regressions, where we used tissue or cell type specific gene expression value as an outcome of a linear regression and binary status of genes as a predictor while conditioning on average expression and gene size). Note that the conditional analyses were performed for gene-sets identified by analysis based on full GWAS results but not for the set of prioritized genes since the enrichment was tested by the hypergeometric test. To systematically identify independent and dependent associations of gene-properties (or gene-sets), and group them into clusters of independent associations, we set thresholds for conditional P-value and proportional significance (PS) of conditional P-value relative to marginal P-value. To simplify, we use cell type A and B for examples hereafter, while same applies to tissue type or gene-set. Here we denote, p_A and p_B as the marginal P-values of cell type A and B where p_A<p_B, p_A,B as the P-value of cell type A conditioning on cell type B, and PS_A as the proportional significance of cell type A conditioning on cell type B (computed as -log10(p_A,B)/-log10(p_A)). We then defined relationship between cell type A and B based on the decision rules as below (ordered according to the priority of the rules).

• When PS_A≥0.8 and PS_B≥0.8, associations are independent.

• When p_A,B≥0.05 and p_B,A≥0.05, both cell types are jointly associated and both are equally likely true signals (grouped into a same cluster).

• When PS_A≥0.5 and p_B,A≥0.05, the association of cell type B is confounded by cell type A and cell type B is discarded.

• When PS_A≥0.5 and PS_B<0.2, the association of cell type B is mostly dependent on cell type A and cell type B is discarded.

• When PS_A≥0.2 and p_B,A≥0.05, the association of cell type A is partially jointly explained by cell type B but the association of cell type B is confounded by A, therefore cell type B is discarded.

• When PS_A≥0.2 and PS_B<0.2, the association of both cell types are partially jointly explained while the association of cell type B is largely depending on cell type A, therefore cell type B is discarded.

• When PS_A≥0.5 and PS_B<0.5, same as above.

• When PS_A<0.2 and PS_B<0.2, while there are still remaining signals by conditioning on each other, both cell types are largely jointly explained (grouped into a cluster).

• When PS_A≥0.2, PS_B≥0.2, p_A,B≥0.01 and p_B,A≥0.01, although there are still small signal remaining in both cell types, associations of both cell types are considered as jointly explained (grouped into a cluster).

• When PS_A≥0.2 and PS_B≥0.2, both cell types are partially jointly explained while there are remaining signals, therefore cell type A and B are grouped into separate clusters.

• When cell type A and B are collinear (in a case when regression fails), they are considered as joint associations and grouped into a cluster.

The conditional analyses were performed with the most significant sets until all cell types are denoted independent, clustered or discarded.