Our individually variable input GWAS data were obtained from 6 GWAS involving aspects related to poor sleep quality, including Trouble falling asleep (TFA), Insomnia (Ins), Undersleep (Und), Sleep disorders (SlD), Hypnotic drug dependence (HDD), and Tiredness (Tir). All input GWAS had ethical clearance from their respective institutional review boards, and all participants provided informed consent, and the data were subjected to strict quality control. Of these, TFA (n = 243,876) was from a study by Verma et al. (2024), Und (n = 110,188) was from a study by Jones et al. (2016), HDD (n = 146,106) from a study by Jiang et al. (2021), Ins (n = 497,539) and SlD (n = 495,270) from Finngen: https://www.Finngen.fi/en (Kurki et al., 2023), Tir (n = 449,019) from the IEU OpenGWAS project: https://gwas.mrcieu.ac.uk (Lyon et al., 2021; detailed GWAS listing information is available in Supplementary Table 1).

Firstly, it is imperative to exclude low-quality samples, defined as those with a missing rate exceeding 5%. Next, the MHC region (MHC, major histocompatibility complex) is located on chromosome 6, at a specific location between approximately 25,000,000 and 35,000,000 base pairs (genomic position). Due to the genetic diversity and structural complexity of the MHC region, especially the polymorphism of immune-related genes (Plasil et al., 2019; Lie and Thorsby, 2005), the MHC region is usually specially processed. Subsequently, the construction of a GWAS was initiated. In this study, the default parameters were utilized during the preparation of summary statistics. The autosomal SNPs from the five input PSQ GWAS that passed the recommended default quality control filters were filtered to the 1,000 Genomes Phase 3 EUR panel. SNPs with a minor allele frequency (MAF) < 0.01 were removed. These SNPs are error-prone due to the small number of samples in the genotype cluster, and the standard error of the regression of the linkage disequilibrium (LD) score for these SNPs is usually high (Lachance, 2010). SNPs with zero effect estimate were removed (to avoid affecting matrix reactivity, which is necessary for the Genomic-SEM), SNPs that did not match the reference panel were removed, palindromic SNPs with uncertain direction were excluded, strand inversions were corrected by referencing a unified allele coding scheme, and SNPs with mismatched alleles were excluded.

In our analysis, the single-input GWAS we included came from different genome repositories and had different participants. This means that when conducting the GWAS, we fully considered the sample overlap between different cohorts to ensure the accuracy and completeness of the results, as well as the statistical impact of potential sample overlap.

We used the GenomicSEM R package (v.0.0.5) to implement genomic structural equation modeling of TFA, Ins, Und, SlD, HDD, and Tir in a GWAS analysis to investigate the broad genetic susceptibility behind these PSQ related traits (Figure 1, By Figdraw). Genomic-SEM is a newly developed multivariate approach that can examine multiple latent multivariate models to explore the potential structure of traits of interest (Grotzinger et al., 2019; refer to Table 1 for detailed criteria.).

Genomic-SEM is not biased by sample overlap (e.g., UKB participants overlapping in multiple input GWAS validation) or sample size imbalance. Furthermore, this method facilitates the identification of variants that affect only some but not all complex traits, which do not represent a broad cross-trait susceptibility.

Genomic-SEM is a two-stage process (Nock and Zhang, 2011). In the first stage, the empirical genetic covariance matrix and the corresponding sampling covariance matrix are estimated. We prepared summary statistics for the PSQ GWAS for the first stage and used a multivariate extension of cross-trait LDSC (Linkage disequilibrium score regression; Bulik-Sullivan et al., 2015) to generate the empirical genetic covariance matrix between the six traits as input for the SEM common factor model. In the subsequent stage, an SEM model was specified with the objective of minimizing the discrepancy between the assumed and the empirically calculated covariance matrices from the initial stage. The primary research objective was to identify the genetic underpinnings of the six sleep quality-related traits. To that end, a univariate model was tested. The model fit was evaluated using SRMR, model χ², the Akaike information criterion, and CFI (comparative fit index). By implementing the appropriate common-factor SEM specification, individual autosomal SNPs were incorporated into the genetic and related sample covariance matrices, thereby yielding a polygenic broad-sense heritability result of 5,761,413 shared covariance between related GWAS.

To assess whether SNP associations are appropriately modeled within a multivariate structural equation modeling (SEM) framework, we computed the SNP heterogeneity statistic. The original hypothesis of the SNP test was that SNP associations in a single phenotypic GWAS are statistically moderated by a constructed GWAS for PSQ. Thus, significant QSNP testing in the constructed GWAS suggests that SNPs may be associated with pathways other than the shared genetic mechanisms in the constructed model. SNPs with significance thresholds (p < 0.05) were subsequently excluded.

The FUMA GWAS (Functional Mapping and Annotation of Genome-Wide Association Studies) method was utilized to identify genomic loci and to ascertain the leading SNP loci associated with the constructed GWAS (Watanabe et al., 2017; Watanabe et al., 2019). These SNPs exhibited a low correlation (less than 0.1) with other SNPs in LD and possessed genome-wide significance (p value < 5 × 10⁻⁸). Initially, the summary statistics of the SNPs from the constructed GWAS were entered into the program to assess their association strength. Furthermore, a comparison was made between the leading SNP loci and the original single-input GWAS. To ascertain whether the 14 leading SNP loci in the new GWAS were associated with multiple effects, the GWAS Catalog was consulted for published significant associations (p value < 5 × 10⁻⁸). Furthermore, a risk gene locus analysis was performed on the established model based on the fuma software function, with a significance threshold of p value < 5 × 10⁻⁸, and the relevant output file was analyzed with MAGMA (Multi-marker Analysis of GenoMic Annotation). MAGMA is a tool for post-processing GWAS that aims to assess associations between genes and phenotypes such as disease or health traits. The MAGMA tool integrates multiple genetic markers (e.g., SNPs) into a gene-level signal, and calculates the association of each gene with the phenotype. The objective is to extract information about gene function from genome-wide SNP data in order to analyze genetic signals at the gene level, with a significance threshold of FDR-p value < 0.05. Furthermore, a novel methodology has been developed, which is referred to as the GWAS locus reduction method. This method compares the genome-wide significance threshold used in single-input GWAS with the lead sites identified by Genomic-SEM. The implementation of this method has the potential to facilitate the identification of additional valuable novel sites at the lead site.

In order to identify the most likely causal variants associated with the new GWAS, SuSIE (Sum of Single Effects) and FINEMAP were utilized, with the latter being implemented in the R package echolocatoR v.2.0.3 (Akdeniz et al., 2024). Identifying possible causal variants using SuSIE and FINEMAP: SuSIE and FINEMAP are both tools for fine-mapping analysis, which aim to identify the most likely causal variants associated with a phenotype. In this step, a 250 kb window was used to include the regions associated with each lead SNP, and the causal inference probability of each SNP within these regions was calculated. Confident set: A probability threshold of 0.95 was set. If the posterior probability of a variant exceeds this threshold, it is designated as a possible causal variant. Consensus SNP and probability set: echolocatoR defines a ‘consensus SNP’, that is, a variant that appears in both SuSIE and FINEMAP results. For these consensus SNPs, the tool calculates their average posterior probability and determines the average credible set based on the probability results. The credibility is defined as 1 when the posterior probability of the SNP in SuSIE and FINEMAP exceeds 0.95, otherwise it is 0.

Following the localisation of possible causal variants, a TWAS (Transcriptome-Wide Association Study) was performed to prioritize genes associated with the constructed GWAS based on the relationship between gene expression and phenotype (Mai et al., 2023; Li et al., 2021). The FUSION method was used for TWAS, and 37,920 pre-computed expression quantitative trait loci (eQTL) traits from GTExv.8 data were utilized. These were then used to calculate the expression associations between different genes and tissues.

A subsequent analysis of the TWAS results revealed that the new GWAS data contained sufficient variation to analyze 36,149 traits (from 37,920 eQTL traits), thereby indicating the data’s high quality. Genes with a p value of less than 0.05 (genes significantly associated with the constructed SEM) were included in further analysis. For these TWAS significant genes, the FOCUS method (a fine-mapping method designed specifically for TWAS studies) was further performed, and the FOCUS method was used to prioritize novel GWAS genes. The FOCUS method assesses whether there is a causal relationship between a gene and a phenotype based on the FOCUS posterior inclusion probability. Combining with previous studies, we considered TWAS significant genes that not only showed significance in the TWAS analysis, but also were consistent with other evidence (such as FOCUS), indicating that they may be causal.

In order to investigate the potential relationship between PSQ and Mendelian disease genes and their related pathways, MAGMA and FUMA (GESA) data were used for gene enrichment analysis and gene pathway set analysis. In addition, MendelVar¹ was used for gene enrichment analysis (Sobczyk et al., 2021).

In order to identify the cell types associated with PSQ, Cell Type Expression Specificity Integration for Single-Cell RNA Sequencing Data of Complex Traits (CELLECT) was utilized (Timshel et al., 2020), with the Tabula Muris90 dataset being employed as a source of data (Tabula Muris Consortium, Overall coordination, Logistical coordination, Organ collection and processing, Library preparation and sequencing, Computational data analysis, Cell type annotation, 2018). This particular dataset contains transcriptome data from 100,000 cells and 20 organs and tissues from mice (Mus musculus). The CELLEX was then utilized for the preprocessing and normalization of the single-cell RNA sequencing data from Tabula Muris, followed by the calculation of the expression-specific likelihood score for each gene. Subsequently, the LDSC software was employed for cell type-specific analysis, and the cell types were classified.

The LDSC tool is utilized to compute partitioning heritability, with the contribution of each genomic region to heritability being assessed by assigning the genetic information of a phenotype to different genomic regions (genes, enhancers, silencers, etc.). Specifically, LDSC employs a weighted LD matrix, a genotype frequency file, and summary statistics to perform the calculation. This process ultimately provides an estimation of the genetic contribution of each region.

In order to identify the degree of association between previously measured diseases and biomarkers and the constructed GWAS data on PSQ that was not directly measured, a large-scale correlation analysis was performed. The analysis incorporated 13,014 phenotypes from the IEU database as potential exposure factors and 20 common phenotypes of neurological diseases from the IEU and FinnGen databases as potential outcome factors. Mendelian randomization (MR) analysis was performed using the TwoSampleMR package v.0.6.8, with inverse variance weighted (IVW) as the primary correlation test (p < 0.05). The final results underwent Bonferroni multiple correction.

Polygenic risk scores (PRS; Ge et al., 2019) were calculated based on genome-wide summary statistics, and the genetic contribution of different chromosomal regions to the development of the disease was assessed. The method utilized PRS-CS (Polygenic Risk Score with Continuous Shrinkage) software to estimate the posterior effect values of SNPs through GWAS data and an external LD reference panel. PRS-CS was implemented using the standard European LD reference panel from the 1,000 Genomes Project, and default PRS-CS parameters were used. PRS is calculated using a Bayesian regression model that integrates an LD reference panel based on GWAS summary statistics.

LD-Score regression analysis of the six univariate inputs into the GWAS for TFA, Ins, Und, SlD, HDD, and Tir had heritability contributions Z-values of 15.9, 5.49, 10.5, 20.1, 1.44, and 23.8, respectively. The Z-values of the genetic covariates for their two bivariate inputs were 10.8 (TFA and Und), 16.2 (TFA and Tir), 2.31 (Ins and HDD), and detailed one-way genetic parameters are detailed in Supplementary Table 2. The six genetic covariance matrices fitted well to the common factor model, as indicated by the following fit indices: comparative fit index (CFI) = 0.93 and standardized root mean square residual (SRMR) = 0.09. A detailed assessment of the model stability is provided in Supplementary Table 3. The univariate SEM parameters for the potential factors (F1) are presented in Supplementary Table 4. These results collectively suggest evidence of shared genetic factors.

By extending SEM to incorporate individual variation, we generated an indirectly measured obtained GWAS estimating the association of 5,761,413 SNPs with PSQ.

We have removed a total of 4,777,855 SNPs through parameter control in the method, and a total of 983,558 effective SNPs have been retained after retaining the regression coefficients. The Mean chi² value for all SNPs is 1.45, the genome control Lambda GC is 1.398, extreme value standard Max Chi² is 40.019, genome-wide significance level is 78, total observed scale heritability (h²) is 0.0038(0.0002), genetic contribution to environmental contribution Ratio is 0.0178(0.0214), intercept term in the regression model is 1.008, standard error of the intercept term in the regression model is 0.0096. The multiple estimates directly indicate that the potential inflation of the structural equation we are concerned with is due to the polygenic heritability signal, rather than population stratification bias, and pleiotropic parameter effects.

Using FUMA software (Figures 2a,b), we evaluated the Genomic-SEM and found 14 leading SNP loci (Figure 2c; Supplementary Table 5a). We also identified a total of 460 potential genes associated with PSQ through genome-wide significant control (Significance was defined as a genome-wide p < 5 × 10⁻⁸; findings with FDR < 0.05 were considered suggestive; Supplementary Table 5b). The vast majority of these 14 leading SNP loci are located in intronic and intergenic (Supplementary Table 6). In addition, a total of 8 GWAS reduction loci (including more valuable novel loci identified among leading SNP loci) were identified (such as rs2820309, rs6421926, rs4588900, etc.; Supplementary Table 7). rs6421926 has been reported in multiple research papers (Akimova et al., 2025; Baselmans et al., 2019), but we found that this site is not directly related to PSQ. Instead, it is a potential mediator site. rs4588900 has been reported in previous studies of European populations to be associated with snoring (Watanabe et al., 2022; Jansen et al., 2019), which is relevant to the phenotypes in our study.

Fine Mapping Analysis identified strong associations at multiple genomic loci (mean. PP > 0.95), including: chromosome 14 (rs8003028, rs2274077, rs77168063, variant in KTN1). Regional plots showed significant peaks at these loci, and other credible set variants also showed evidence of association (Figure 3a; Supplementary Table 8).

Next, we performed a Transcriptome-Wide Association Study (TWAS) using FUSION to identify gene-level associations with PSQ. We found 23 genes that passed multiple comparison correction (Supplementary Table 9; Figure 3b). Next, we performed a fine mapping analysis using FOCUS on genomic structural equation data, and 18 genes were found to be possible disease-causing signals for PSQ (pips>0.8; Supplementary Table 10). To further confirm these “highly credible” gene-level associations, we performed an intersection test. These include LMOD1, ZDHHC5, MAD1L1, MED19, etc. (Figure 3c). Among them, the TWAS Z-scores of LMOD1, ZDHHC5, etc., were all greater than 0, indicating that the predicted gene expression was positively correlated with PSQ, suggesting that upregulation of these genes may be associated with PSQ. In contrast, the TWAS Z-scores of MAD1L1, MED19, etc., were less than 0, indicating that downregulation of these genes may be associated with PSQ (Supplementary Table 11; Figure 4a).

Multi-marker Analysis of GenoMic Annotation (MAGMA) identified 35 genes (bon-p < 0.05; Supplementary Table 12; Figure 4b), which we utilized for gene set analyses, and these genes showed enrichment in GSEA entries (Supplementary Table 13); many of the gene sets were associated with Glutamate synapse, γ-Aminobutyric acid synapse, Neuroticism, etc. In addition, biological processes mapped by MendelVar enrichment were supported by mapping in GSEA entries (autosomal dominant intellectual developmental disorder 4), however, there were no significant results after FDR correction (Figure 5). Analysis of enrichment from different cell types showed 5 cell types exceeding the significance criterion (Supplementary Table 14). Four of the 5 cell types were associated with the nervous system, namely Neuron, Oligodendrocyte, Oligodendrocyte precursor cell and Astrocyte. After FDR correction Neuron still exceeded the significance criterion (FDR-p < 0.05).

In the results of the contribution of heritability on genomic regions, we found that most of the genetically contributing sites are concentrated in the regulatory regions of chromosomes and histone modification regions such as H3K4me1. These regions are usually key sites for gene expression regulation, chromatin modification, and transcription factor binding. In particular, the effects of genetic variation are most pronounced in regions of conserved genomic regions and intronic regions, which may play an important role in traits or disease susceptibility by regulating gene expression levels. In addition, certain non-coding regions such as the 3’ Untranslated Region also showed strong genetic contributions, suggesting that these regions may be involved in complex genetic mechanisms by regulating gene expression or function (Supplementary Table 15).

In the results of biomarker analysis we found 705 positive exposures (Supplementary Table 16), in which genes such as PDK4 (ENSG00000004799) were positively correlated with PSQ, and on the contrary genes such as MAD1L1 (ENSG00000002822) were negatively correlated with it, and traits such as Daytime nap, Body mass index (BMI) may be risk factors for PSQ, while such as HDL cholesterol levels may be protective factors. After Bonferroni multiple correction, there are still 93 positive exposure factors. All of these factors may have a biological impact on sleep quality. In addition, we further investigated the potential relationship between PSQ and 20 common neurological disorders and found that PSQ was a risk factor for Stroke and Ischemic stroke, and all of them passed the tests of heterogeneity and multiplicity (Supplementary Table 17; Figure 6).

Our analyses showed that PRS preformed variant loci were strongly associated with the risk of developing the disease and that the genetic contribution to the disease varied significantly among different chromosomal regions. In particular, we observed higher genetic contributions in regions such as chromosome 6 (−0.061) and chromosome 8 (−0.063), which may contain important genes and regulatory elements that influence disease susceptibility (Supplementary Table 18).