Ethics approval
All of Us (All of Us Institutional Review Board), UK Biobank (North West Multi-centre Research Ethics Committee/North West–Haydock Research Ethics Committee), FinnGen (Coordinating Ethics Committee of the Hospital District of Helsinki and Uusimaa), Estonian Biobank (Estonian Committee on Bioethics and Human Research, Estonian Ministry of Social Affairs), Genes & Health (London–South East Research Ethics Committee/NRES Committee London–South East), Mass General Brigham Biobank (Mass General Brigham Institutional Review Board/Mass General Brigham Human Research Committee), Michigan Genomics Initiative (University of Michigan Medical School Institutional Review Board), deCODE Genetics/Amgen Iceland (National Bioethics Committee of Iceland and Icelandic Data Protection Authority), Copenhagen Hospital Biobank (CHB; Danish National Committee on Health Research Ethics and Capital Region Data Protection Agency), Danish Blood Donor Study (DBDS; Danish National Committee on Health Research Ethics and Danish Data Protection Agency), Intermountain Health (Intermountain Healthcare Institutional Review Board), and Nashville Biosciences/BioVU (Vanderbilt University Medical Center Institutional Review Board) gave ethical approval for this work. Participants were not compensated for this study.
Phenotype definition
Cases were individuals with at least one inpatient or primary care record containing ICD-10 code M79.7 (fibromyalgia; Supplementary Table 11). No exclusion criteria were applied: controls were all genotyped participants without code M79.7. Phenotype ascertainment was performed independently within each cohort before genome-wide association analysis.
Genotyping, imputation and quality control
We included autosomal and chromosome X data for all cohorts.
All of Us
We obtained whole-genome sequencing (WGS), hg38 data from All of Us’s Curated Data Repository version 8 release, which involved variant calling, initial quality control and genetic ancestry assignment. We used the allele count/allele frequency (ACAF) threshold callset, which contains variants with a minor allele count >100 or minor allele frequency >1% in any genetic ancestry group, as the basis for a second round of ancestry-specific quality control with version 2.0.0 of the plink genetic analysis toolkit62, in which variants with minor allele frequency (MAF) < 0.1%, missingness >10% or Hardy–Weinberg equilibrium (HWE) P < 1 × 10−15 (with mid-P correction) were removed. Participants flagged by All of Us in the first round of quality control, those with sex-chromosome aneuploidy or sex at birth not equal to male or female, were excluded. We performed association analysis on each of the 3 largest genetic ancestry groups—EUR, AFR and AMR—with version 3.3 of the regenie GWAS toolkit63 applying Firth approximation for variants with association value < 0.01, and using as covariates age, sex, age2, age × sex, age2 × sex and the first 10 genotype principal components (PCs). For step 1 of regenie, we used an LD-pruned subset of the full genotypes, calculated with plink version 2.0.0 via the option ‘–indep-pairwise 500 kb 1 0.2’.
We gratefully acknowledge All of Us participants for their contributions, without whom this research would not have been possible. We also thank the National Institutes of Health’s All of Us Research Program for making available the participant data examined in this study.
UK Biobank
We obtained imputed genotypes from the UK Biobank’s Data-Field 22828. The UK Biobank imputed genotypes using a combination of two GRCh37 reference panels, the Haplotype Reference Consortium and a combined UK10K and 1000 Genomes Phase 3 panel, using the Haplotype Reference Consortium imputation if a variant was present in both panels64. We excluded samples with sex-chromosome aneuploidy (Data-Field 22019), that were outliers for heterozygosity or missingness rate (Data-Field 22027), or that had discordant genetic sex (Data-Field 22001) versus self-reported sex (Data-Field 31). We excluded variants with MAF < 0.1%, imputation INFO score <0.8, missingness >10% or HWE P < 1 × 10−15 (with mid-P correction) with plink version 2.0.0. We performed associations on the 3 largest genetic ancestry groups as defined by Pan-UKBB (Return 2442)—EUR, Central and South Asian (CSA) and AFR—with version 3.4.1 of regenie, applying Firth approximation for variants with association P value < 0.01 and using as covariates age, sex, age2, age × sex, age2 × sex and the first 10 genotype PCs. For step 1 of regenie, we used an LD-pruned subset of the full genotypes, calculated with the ‘–indep-pairwise 500 kb 1 0.2’ option from plink version 2.0.0.
We thank all participants of the UK Biobank. This research has been conducted using the UK Biobank Resource under application number 116200.
FinnGen
Genotyping in the FinnGen cohort was performed by using Illumina and Affymetrix arrays (Thermo Fisher Scientific) and lifted over to GRCh38/hg38. Individuals with high genotype absence (>5%), inexplicit sex or excess heterozygosity (±4 standard deviations) were excluded from the data. In addition, variants that had high absence (>2%), low minor allele count (<3) or low HWE (P < 1 × 10−6) were removed. More detailed explanations of the genotyping, quality control and genotype imputation are provided elsewhere65. All individuals in the cohort were Finns and matched against the SiSu v4 reference panel.
For the FinnGen cohort (Data Freeze 12), GWAS was conducted using the regenie pipeline (https://github.com/FINNGEN/regenie-pipelines). Analysis was adjusted for age at death or end of follow-up, sex, genotyping batches and the first ten genetic PCs. Firth approximation was applied for variants with association P value < 0.01.
Study subjects in FinnGen provided informed consent for biobank research, based on the Finnish Biobank Act. Alternatively, separate research cohorts, collected before the Finnish Biobank Act came into effect (in September 2013) and start of FinnGen (August 2017), were collected based on study-specific consents and later transferred to the Finnish biobanks after approval by the Finnish Medicines Agency and the National Supervisory Authority for Welfare and Health. Recruitment protocols followed the biobank protocols approved by the Finnish Medicines Agency. The Coordinating Ethics Committee of the Hospital District of Helsinki and Uusimaa (HUS) statement number for the FinnGen study is Nr HUS/990/2017.
The FinnGen study is approved by the Finnish Institute for Health and Welfare (permit numbers THL/2031/6.02.00/2017, THL/1101/5.05.00/2017, THL/341/6.02.00/2018, THL/2222/6.02.00/2018, THL/283/6.02.00/2019, THL/1721/5.05.00/2019 and THL/1524/5.05.00/2020), Digital and Population Data Service Agency (permit numbers VRK43431/2017-3, VRK/6909/2018-3, VRK/4415/2019-3), the Social Insurance Institution (permit numbers KELA 58/522/2017, KELA 131/522/2018, KELA 70/522/2019, KELA 98/522/2019, KELA 134/522/2019, KELA 138/522/2019, KELA 2/522/2020 and KELA 16/522/2020), Findata (permit numbers THL/2364/14.02/2020, THL/4055/14.06.00/2020, THL/3433/14.06.00/2020, THL/4432/14.06/2020, THL/5189/14.06/2020, THL/5894/14.06.00/2020, THL/6619/14.06.00/2020, THL/209/14.06.00/2021, THL/688/14.06.00/2021, THL/1284/14.06.00/2021, THL/1965/14.06.00/2021, THL/5546/14.02.00/2020, THL/2658/14.06.00/2021 and THL/4235/14.06.00/2021), Statistics Finland (permit numbers TK-53-1041-17 and TK/143/07.03.00/2020 (earlier TK-53-90-20) TK/1735/07.03.00/2021, TK/3112/07.03.00/2021) and Finnish Registry for Kidney Diseases permission and extract from the meeting minutes on 4 July 2019.
The Biobank Access Decisions for FinnGen samples and data used in FinnGen Data Freeze 12 include THL Biobank BB2017_55, BB2017_111, BB2018_19, BB_2018_34, BB_2018_67, BB2018_71, BB2019_7, BB2019_8, BB2019_26, BB2020_1 and BB2021_65; Finnish Red Cross Blood Service Biobank 7.12.2017; Helsinki Biobank HUS/359/2017, HUS/248/2020, HUS/430/2021 §28 and §29, HUS/150/2022 §12, §13, §14, §15, §16, §17, §18, §23, §58 and §59, and HUS/128/2023 §18; Auria Biobank AB17-5154 and amendment 1 (17 August 2020) and amendments BB_2021-0140, BB_2021-0156 (26 August 2021, 2 Feb 2022), BB_2021-0169, BB_2021-0179, BB_2021-0161, AB20-5926 and amendment 1 (23 April 2020) and its modifications (22 September 2021), BB_2022-0262, BB_2022-0256, Biobank Borealis of Northern Finland_2017_1013, 2021_5010, 2021_5010 Amendment, 2021_5018, 2021_5018 Amendment, 2021_5015, 2021_5015 Amendment, 2021_5015 Amendment_2, 2021_5023, 2021_5023 Amendment, 2021_5023 Amendment_2, 2021_5017, 2021_5017 Amendment, 2022_6001, 2022_6001 Amendment, 2022_6006 Amendment, 2022_6006 Amendment, 2022_6006 Amendment_2, BB22-0067, 2022_0262, 2022_0262 Amendment, Biobank of Eastern Finland 1186/2018 and amendment 22§/2020, 53§/2021, 13§/2022, 14§/2022, 15§/2022, 27§/2022, 28§/2022, 29§/2022, 33§/2022, 35§/2022, 36§/2022, 37§/2022, 39§/2022, 7§/2023, 32§/2023, 33§/2023, 34§/2023, 35§/2023, 36§/2023, 37§/2023, 38§/2023, 39§/2023, 40§/2023 and 41§/2023; Finnish Clinical Biobank Tampere MH0004 and amendments (21.02.2020 and 06.10.2020), BB2021-0140 8§/2021, 9§/2021, §9/2022, §10/2022, §12/2022, 13§/2022, §20/2022, §21/2022, §22/2022, §23/2022, 28§/2022, 29§/2022, 30§/2022, 31§/2022, 32§/2022, 38§/2022, 40§/2022, 42§/2022 and 1§/2023; Central Finland Biobank 1-2017, BB_2021-0161, BB_2021-0169, BB_2021-0179, BB_2021-0170, BB_2022-0256 and BB_2022-0262, BB22-0067; decision allowing to continue data processing until 31 August 2024 for projects BB_2021-0179, BB22-0067,BB_2022-0262, BB_2021-0170, BB_2021-0164, BB_2021-0161 and BB_2021-0169; Terveystalo Biobank STB 2018001 and amendment 25 August 2020; Finnish Hematological Registry and Clinical Biobank decision 18 June 2021; and Arctic biobank P0844: ARC_2021_1001.
Estonian Biobank
All the Estonian Biobank participants have been genotyped at the Core Genotyping Lab of the Institute of Genomics, University of Tartu, using Illumina Global Screening Array (versions 1, 2 or 3). Samples were genotyped and PLINK format files were created using Illumina GenomeStudio v2.0.4. Individuals were excluded from the analysis if their call rate was less than 95%, if they were outliers of the absolute value of heterozygosity (>3 standard deviations from the mean) or if sex defined based on heterozygosity of the X chromosome did not match sex in phenotype data. Before imputation, variants were filtered by call rate <95%, HWE P < 1 × 10−4 (autosomal variants only) and minor allele frequency <1%. Genotyped variant positions were in build 37 and were lifted over to build 38 using Picard. Phasing was performed using the Beagle v5.4 software. Imputation was performed with Beagle v5.4 software (22 July 2022 release) and default settings. The dataset was split into batches of 5,000. A population-specific reference panel consisting of 2,695 WGS samples was used for imputation, and standard Beagle hg38 recombination maps were used. On the basis of the PC analysis, samples that were not of European ancestry were removed. Duplicate and monozygous twin detection was performed with KING 2.2.7, and one sample was removed from the pair of duplicates. Analyses were restricted to individuals of European ancestry.
Association analysis in the Estonian Biobank was carried out for all variants with an INFO score >0.4 using the additive model as implemented in regenie v3.0.3 with standard binary trait settings. Logistic regression was carried out with adjustment for current age, age2, sex and the first 10 genetic PCs as covariates, analyzing only variants with a minimum minor allele count of 2.
The activities of the Estonian Biobank are regulated by the Human Genes Research Act, which was adopted in 2000 specifically for the operations of the Estonian Biobank. Individual-level data analysis in the Estonian Biobank was carried out under ethical approval 1.1-12/624 from the Estonian Committee on Bioethics and Human Research (Estonian Ministry of Social Affairs), using data according to release application 3-10/GI/31689 from the Estonian Biobank.
We want to acknowledge the participants of the Estonian Biobank for their contributions. The Estonian Genome Center analyses were partially carried out in the High Performance Computing Center, University of Tartu. The Estonian Biobank Research Team was responsible for data collection, genotyping, quality control and imputation and consisted of A. Metspalu (andres.metspalu@ut.ee), M. Metspalu (mait.metspalu@ut.ee), L. Milani (lili.milani@ut.ee), R. Mägi (reedik.magi@ut.ee), M. Nelis (mari.nelis@ut.ee), T. Esko (tonu.esko@ut.ee) and G. Hudjashov (georgi.hudjashov@ut.ee).
Genes & Health
Genotyping was performed on Illumina Infinium Global Screening Array v3 with additional multi-disease variants. Variants with call rates less than 0.99 and/or MAF < 1% were excluded. We excluded individuals unlikely to have genetically inferred Pakistani or Bangladeshi ancestry. Imputation was performed using the TOPMed-r2 panel. We excluded single nucleotide polymorphisms (SNPs) with low imputation scores (INFO < 0.3). In the Genes & Health cohort, sex was defined on the basis of XX (female) and XY (male) chromosomal presence in genotype data.
Diagnoses were curated from routine UK NHS EHR data from primary care (Systematized Nomenclature of Medicine coded) and secondary care (ICD-10 coded) sources. Data were combined without mapping between coding formats. For each clinical code, the earliest ever measure recorded in a participant’s medical records was used, excluding erroneous code dates preceding the participant’s recorded date of birth. Association analysis was run in regenie, applying Firth approximation for variants with association P value < 0.01 and using as covariates age, sex, age2 and the first 20 genotype PCs.
Genes & Health is and has recently been core funded by Wellcome (WT102627, WT210561), the Medical Research Council (UK) (M009017, MR/X009777/1, MR/X009920/1), Higher Education Funding Council for England Catalyst, Barts Charity (845/1796), Health Data Research UK (for London substantive site) and research delivery support from the NHS National Institute for Health Research Clinical Research Network (North Thames). We acknowledge the support of the National Institute for Health and Care Research Barts Biomedical Research Centre (NIHR203330); a delivery partnership of Barts Health NHS Trust, Queen Mary University of London, St George’s University Hospitals NHS Foundation Trust and St George’s University of London.
Genes & Health is and has recently been funded by Alnylam Pharmaceuticals, Genomics and a Life Sciences Industry Consortium of AstraZeneca, Bristol-Myers Squibb, GlaxoSmithKline Research and Development, Maze Therapeutics, Merck Sharp & Dohme, Novo Nordisk A/S, Pfizer and Takeda Development Centre Americas.
We thank Social Action for Health, Centre of The Cell, members of our Community Advisory Group and staff who have recruited and collected data from volunteers. We thank the NIHR National Biosample Centre (UK Biocentre), the Social Genetic & Developmental Psychiatry Centre (King’s College London), Wellcome Sanger Institute and Broad Institute for sample processing, genotyping, sequencing and variant annotation. This study uses data provided by patients and collected by the NHS as part of their care and support. This research used Queen Mary University of London’s Apocrita HPC facility, supported by QMUL Research-IT, https://doi.org/10.5281/zenodo.438045.
We thank Barts Health NHS Trust, NHS Clinical Commissioning Groups (City and Hackney, Waltham Forest, Tower Hamlets, Newham, Redbridge, Havering, Barking and Dagenham), East London NHS Foundation Trust, Bradford Teaching Hospitals NHS Foundation Trust, Public Health England (especially D. Wyllie), Discovery Data Service and Endeavour Health Charitable Trust (especially D. Stables), Voror Health Technologies (especially S. Don), NHS England (for what was NHS Digital) for GDPR-compliant data sharing backed by individual written informed consent.
Most of all, we thank all of the volunteers participating in Genes & Health.
A favorable ethical opinion for the main Genes & Health research study was granted by NRES Committee London – South East (reference 14/LO/1240) on 16 September 2014. Queen Mary University of London is the sponsor and data controller.
Mass General Brigham Biobank
The Mass General Brigham Biobank genotyped 53,297 participants on the Illumina Global Screening Array (GSA) and 11,864 on Illumina Multi-Ethnic Global Array (MEG). The GSA arrays captured approximately 652,000 SNPs and short insertions–deletions (indels), while the MEG arrays captured approximately 1.38 M SNPs and short indels. These genotypes were filtered for high missingness (>2%) and variants out of HWE (P < 1 × 10−12), as well as variants with an allele frequency (AF) discordant (P < 1 × 10−150) from a synthesized AF calculated from gnomAD subpopulation frequencies and a genome-wide GnomAD model fit of the entire cohort. This resulted in approximately 620,000 variants for the GSA and 1.15 M for MEG. The two sets of genotypes were then separately phased and imputed on the TOPMed imputation server (Minimac4 algorithm) using the TOPMed-r2 reference panel. The resultant imputation sets were both filtered at an r2 > 0.4 and an MAF > 0.001, and then the two sets were merged and intersected resulting in approximately 19.5 M hg38 autosomal variants. The sample set for analysis was then restricted to just those classified as EUR according to a random-forest classifier trained with the Human Genome Diversity Project66 as the reference panel, with the minimum probability for assignment to an ancestral group of 0.5, in 19 of 20 iterations of the model. To correct for population stratification, PCs were computed in genetically European participants. Association analysis for the full cohort (N = 51,053) was performed with variants using regenie (v3.2.8)63 with adjustment for age, sex, genotype-chip, tranche and the first ten genetic PCs. The summary statistics were filtered with a minimum minor allele count of 50. These summary statistics were then lifted over to human genome version 19 (hg19).
In the MGB Biobank, sex-specific association analyses were performed in the hg38 build and carried out using regenie v3.2.8 with covariates of age, tranche, genotype-chip and the first 10 PCs of ancestry calculated by performing the within-EUR principal component analysis (PCA). The summary statistics were filtered with a minimum minor allele count of 25.
Michigan Genomics Initiative
We included autosomal and chromosome X data from the Michigan Genomics Initiative (MGI), a health system-based biobank of patients recruited primarily during surgical encounters at Michigan Medicine. As of Freeze 6, MGI comprised 80,529 participants with linked genotype and electronic health record data. Detailed descriptions of the MGI cohort, recruitment protocols and overall design are available in a previous study67.
DNA samples were genotyped at the University of Michigan Advanced Genomics Core on customized versions of the Illumina Infinium CoreExome-24 (v1.0, v1.1, v1.3; ~570,000 markers) or Illumina Infinium GSA (v1.3; ~682,000 markers). Array content included standard backbones with additional custom content to capture GWAS candidate variants, predicted loss-of-function alleles, ancestry-informative markers and pharmacogenomic variants. Genotype calling was performed in GenomeStudio, supplemented by zCall for recovery of rare variants.
Sample-level quality control (QC) excluded individuals for consent withdrawal, genotype-inferred sex mismatch, sex-chromosome aneuploidy, call rate <99%, contamination >2.5%, unresolved technical duplicates or batch-level DNA extraction issues. Relatedness was estimated with KING v2.1.3, and contamination with VICES. Variant-level QC removed probes that did not uniquely map to hg38, sites with call rate <98%, HWE P < 1 × 10−4 in unrelated European-ancestry samples or poor clustering metrics (GenTrain <0.15, cluster separation <0.3). Additional harmonization excluded variants with large allele frequency deviations compared with 1000 Genomes reference populations.
Phasing was performed with Eagle v2.4 using the TOPMed reference panel, followed by imputation against TOPMed haplotypes via the Michigan Imputation Server. Post-imputation, variants with r2 < 0.3 or MAF < 0.01% were removed, yielding ~52 million high-quality variants. PCs were calculated using FlashPCA2 after pruning rare and correlated variants.
For GWAS, we analyzed participants of genetically inferred European ancestry as defined by PCA projection and ADMIXTURE (K = 7 reference populations from the Human Genome Diversity Project). Association testing was performed with SAIGE v0.35, a generalized mixed-model framework that accounts for relatedness and case–control imbalance. Covariates included age, sex, genotyping array and the first ten genetic PCs. Variants with MAF ≥ 0.01% and imputation r2 ≥ 0.3 were included in association testing.
We acknowledge the Michigan Genomics Initiative participants, AI & Digital Health Innovation at the University of Michigan, the University of Michigan Medical School Central Biorepository and the University of Michigan Advanced Genomics Core for providing data and specimen storage, management, processing and distribution services, and the Center for Statistical Genetics in the Department of Biostatistics at the School of Public Health for genotype data curation, imputation and management in support of the research reported in this Article.
deCODE Genetics/Amgen Iceland
At the time of analysis, 63,118 samples from Icelandic participants were whole-genome sequenced at deCODE using Illumina standard TruSeq methods to a mean depth of 38× (ref. 68). Only samples with a genome-wide average coverage of 20× and higher were included. Genotypes of SNPs and indels were identified and called jointly by Graphtyper (v.2.7.1)69. In all, 173,025 samples from Icelandic participants had been chip-genotyped using various Illumina SNP arrays68. The chip-typed individuals were long-range phased70, and the variants identified in the whole-genome sequences of Icelanders imputed into the chip-typed individuals. Using extensive and encrypted Icelandic genealogy data, familial imputation of genotypes in first- and second-degree relatives was used to increase sample size68,70. The final dataset used included 19,657,761 variants with imputation information over 0.8 and MAF over 0.1%.
The Icelandic samples and diagnostic data were obtained from Icelandic medical record data repositories and analyzed under approval from the National Bioethics Committee (17-035-V11-S2, previously 12-162) following review by the Icelandic Data Protection Authority. Data were anonymized and encrypted by a third-party system, approved and monitored by the Icelandic Data Protection Authority.
CHB and the DBDS study
We obtained genotypes from the CHB71 study on pain and degenerative musculoskeletal diseases (CHB-PDS; approval: NVK-1803812, P-2019-51) and the DBDS (approval: NVK-1700407, P-2019-99)72. All samples were genotyped on Illumina’s Infinium Global Screening Array (versions 1.0 and 3.0) and imputed to build hg38. Genotyping and imputation were performed by deCODE Genetics. Before imputation, duplicate samples and those with genotype call rates <98% were removed. Phasing was carried out with SHAPEIT4 (ref. 73), and imputation was performed using deCODE’s in-house workflow, with a joint Graphtyper-based reference panel comprising ~50,000 individuals, including ~10,800 Danes, including danmac5.dk (ref. 74). Initial quality control of genotypes included removal of samples with sex discrepancies, missingness >5% and variants with missingness >10%, or HWE P < 1 × 10−5. For genome-wide association analyses, we removed variants with MAF < 0.1%, imputation INFO < 0.8, missingness >10% or HWE P < 1 × 10−15 (with mid-P correction), using plink version 2.0.0.
Intermountain Health
Under research collaboration between Intermountain Health and deCODE Genetics/Amgen, samples were obtained from consenting participants of European descent in two ongoing studies: the Intermountain Inspire Registry and the HerediGene Population Study75. The Intermountain Healthcare Institutional Review Board approved both studies, and all participants provided written informed consent before enrollment. Samples were genotyped at deCODE Genetics/Amgen using Illumina Global Screening Array chips. In all, 138,006 individuals of European origin were chip typed. Intermountain and Danish imputation was based on a multi-ethnic reference panel of 50,179 whole-genome-sequenced individuals of mostly European descent, including 23,288 individuals from Intermountain and 11,722 individuals from Denmark.
Samples were filtered on 98% variant yield and duplicates removed. The WGS protocol was the same as described above for the Icelandic and Danish data. Sequence variants were imputed into 138,006 chip-typed individuals. Over 245 million high-quality sequence variants and indels, sequenced to a mean depth of 20×, were identified using Graphtyper (v.2.7.1)69. Quality-controlled chip genotype data were phased using SHAPEIT4 (ref. 73). A phased haplotype reference panel was prepared from the sequence variants using the long-range phased chip-genotyped samples68. In all, 21,316,504 variants (imputation info > 0.8 and MAF > 0.1%) were tested.
Nashville Biosciences
Nashville Biosciences (NashBio) is a data and analytics provider owned by Vanderbilt University Medical Center (VUMC). NashBio uses the biobank collection BioVU of VUMC that includes a collection of de-identified DNA samples linked to de-identified data from the electronic health records of VUMC referred to as the Synthetic Derivative (SD) database. All patients of VUMC consented to their residual samples from routine clinical testing and data being contributed to BioVU.
The NashBio dataset included in this study is a subset of BioVU, consisting of, in total, 80,965 whole-genome-sequenced individuals of genetically defined European descent and 31,025 individuals of genetically defined African descent. All samples were whole genome sequenced at deCODE Genetics/Amgen using Illumina NovaSeq in accordance with deCODE/Amgen protocols (deCODE Genetics/Amgen Iceland). Using the same quality criteria as for other analyzed deCODE/Amgen datasets, in total, 22,316,589 variants were tested in the European cohort and 33,108,534 in the African cohort.
BioVU extracts and banks germline DNA samples that are de-identified and only linked to the SD through a randomly assigned unique identifier, not back to the patient or their underlying medical record. The use of BioVU is classified as non-human subject research by the Institutional Review Board (IRB) of VUMC, and NashBio is not required to seek study-specific consent for use of these datasets. The overall biobanking program is reviewed annually by the IRB to maintain this determination and make decisions about patient protections, privacy and ethical issues. Each individual study seeking to use the SD database and BioVU biobank is filed with the IRB of VUMC to validate its non-human subject classification and appropriate use of data.
Genetic ancestry analysis in deCODE datasets
For the non-Icelandic datasets analyzed at deCODE Genetics/Amgen (Copenhagen Hospital Biobank and the DBDS Intermountain Health and Nashville Biosciences), genetic ancestry analysis was performed to identify a subset of individuals with similar ancestry. For the Danish samples, we used ADMIXTURE (v1.23)76 run in supervised mode using the 1000 Genomes populations CEU (Utah residents with northern and western European ancestry), CHB (Han Chinese in Beijing, China), ITU (Indian Telugu in the UK), PEL (Peruvian in Lima, Peru) and YRI (Yoruba in Ibadan, Nigeria) as training samples. These training samples had themselves been filtered for ancestry outliers using PCA and unsupervised ADMIXTURE. Samples assigned <0.93 CEU were excluded, resulting in 358,483 samples included in the analysis. For the Intermountain samples, we used the same methods as for the Danish samples to identify a subset of 108,149 individuals of European descent that were included in the analysis.
To identify ancestry-homogeneous subsets in the Nashville Biosciences cohort, sample ancestry was estimated using similar methods as in the Danish cohort. Samples with admixed African and European ancestry (defined as YRI + CEU > 0.9, YRI > 0.3 and CEU > 0.02) or predominantly African ancestry (YRI > 0.9) were placed in the AFR subcohort. Genotypes were prepared for PCA with PLINK v1.9 (ref. 62) using–maf 0.01–thin-count 1000000–indep-pairwise 60000 6000 0.3 to minimize effects of local LD and very recent population structure. Pairwise relatedness was calculated using KING v2.3.0–ibdseg–degree 3 (ref. 77) to identify relatives of third degree or closer, with one sample from each pair of relatives removed before PCA and projected onto the resulting PCs using a deCODE in-house script that adjusts for shrinkage. PCA was performed with PCAone v0.3.4 (ref. 78) using default parameters.
Association testing in deCODE datasets
The four fibromyalgia case–control GWASs were performed at deCODE Genetics/Amgen (Iceland, Denmark, US Intermountain Health and US NashBio) with software developed at deCODE Genetics, using logistic regression assuming an additive model68. For the Icelandic data, the model included sex, county of birth, current age or age at death (first- and second-order terms included), blood sample availability for the individual, sequencing status and an indicator function for the overlap of the lifetime of the individual with the time span of phenotype collection. To include imputed but ungenotyped individuals in Iceland, we used county of birth as a proxy covariate for the first PC because county of birth has been shown to be in concordance with the first PC in Iceland79. For the Danish data, 12 PCs were used, in addition to sex, year of birth and sequencing status as covariates, whereas for the Intermountain Health data, we included 4 PCs, year of birth, sex and sequencing status as covariates. For the NashBio data, associations were adjusted for the top 20 PCs in addition to sex, year of birth and sequencing batch. The number of PCs used to adjust for population stratification was determined by identifying the point at which further PCs appeared to capture local LD rather than population structure as reflected by sharp peaks in PC loadings and in plateauing of eigenvalues80. All statistical tests were two sided unless otherwise indicated.
Harmonization and meta-analysis
Each cohort’s summary statistics were harmonized via an in-house analysis pipeline. For each summary statistics file (that is, for each cohort and ancestry), we removed variants with missing data in any column, non-ACGT alleles, P values outside (0, 1], allele frequencies outside (0, 1), imputation INFO scores outside (0, 1], or non-positive odds ratios, standard errors, or sample sizes. We normalized indels to their minimal representation by iteratively removing common trailing then leading base pairs from the reference and alternate alleles—subject to the constraint that each allele retains at least one base pair—and advancing the genomic position by the number of leading base pairs removed. We inferred Reference SNP cluster ID (rs) numbers based on chromosome and base-pair columns using dbSNP build 156, requiring an exact match to both the ref and the alt allele, allowing allele flips for single-nucleotide variants (but not indels) and converting the dbSNP variants to minimal representation as well before doing the matching. We removed variants without rs numbers in dbSNP, or that mapped to multiple genomic positions (this happens very rarely). To account for the tendency of dbSNP to merge rs numbers over time, if a variant mapped to multiple rs numbers, we took the lowest-numbered rs number. We calculated each cohort and ancestry’s effective sample size per variant via the formula Neff = ((4/(2 × AAF × (1 − AAF) × INFO)) − BETA2)/SE2, skipping the ‘× INFO’ if INFO was not available, where INFO denotes the imputation information score, BETA the estimated log odds ratio for the variant, and SE its standard error. We then summed these Neff across cohorts and ancestries.
To account for P value inflation due to residual confounding such as cryptic relatedness, we applied stratified LD score regression81 to each cohort and ancestry’s summary statistics before meta-analysis, using the effective sample size Neff rather than the raw N to avoid bias82. We used the standard 52-annotation baseline model precomputed by the developers of stratified LD score regression. We manually calculated LD scores based on ancestry-specific reference panels from the corresponding 1000 Genomes Phase 3 superpopulation, where 1000 Genomes EUR was used for European cohorts, SAS (South Asian) for Central/South Asian cohorts, AFR for African cohorts, and AMR for admixed American cohorts. However, we subset to variants present in Europeans in HapMap 3, as the standard 52 annotations are available only for those variants. For cohort–ancestry pairs in which the LD score regression intercept was greater than 1, we divided each variant’s χ2 statistic by the LD score regression intercept, then adjusted the standard errors and P values accordingly. This has the effect of reducing the significance of every variant in summary statistics with evidence of inflation, while leaving non-inflated summary statistics alone.
After harmonization and LD score regression correction, we performed standard fixed-effects inverse-variance weighted meta-analyses via the ‘–meta-analysis’ option from plink version 1.9.0. Variants were matched across cohorts based on the combination of rs number, reference allele and alternate allele. In total, 54,629 cases and 2,509,126 controls went into the multi-ancestry meta-analysis, and 49,000 cases and 2,264,287 controls went into the European genetic ancestry-only analysis. We also performed sex-stratified meta-analyses restricted to males and females, and leave-one-cohort-out analyses for each of the 11 constituent cohorts, with both multi-ancestry (leaving out all ancestries from that cohort, if the cohort had multiple) and European-only versions. For females, we performed both European-only and multi-ancestry meta-analyses, but for males, we performed only the European-only meta-analysis owing to the lack of male summary statistics in Genes & Health and the low number of male cases in non-European ancestries in the multi-ancestry cohorts.
After meta-analysis, we filtered to variants with an overall minor allele frequency of >1% across all cohorts that went into that particular meta-analysis. As two cohorts (the UK Biobank and Mass General Brigham Biobank) reported variants in hg19 coordinates, we avoided having to perform an (error-prone) remapping of variant positions from hg19 to hg38 by filtering to variants present in at least 5 cohorts in the multi-ancestry meta-analyses and at least 3 studies in the European meta-analyses, thereby ensuring that every variant would be present in at least one hg38 study by the pigeonhole principle.
Causal gene prioritization
We prioritized causal genes for each risk locus in the multi-ancestry meta-analysis via a combination of approaches. First, we performed manual literature search for each lead variant, focusing heavily on the nearest gene to each lead variant as nearest genes are expected to be causal about two-thirds of the time83. We supplemented this search with two orthogonal bioinformatic approaches: the state-of-the-art machine learning method FLAMES84 and an in-house quantitative trait locus (QTL)-based pipeline developed by deCODE Genetics.
FLAMES
FLAMES is a framework that nominates candidate causal genes at risk loci from a GWAS, via an ensemble of two complementary approaches.
First, a gradient boosting model scores genes based on locus-specific variant-to-gene evidence, aggregating features such as QTLs, variant effect predictor annotations, chromatin interactions and variant-to-gene distance. The model’s predictive weights were established by training it on external GWAS loci, in which ‘silver-standard’ causal genes had been identified via missense or predicted loss-of-function variants from exome studies. For this analysis, we supplied the model with 99% credible sets derived from the sum of single effects85 fine-mapping of variants within 500 kb of the lead variant at each locus from our primary meta-analysis.
Second, FLAMES incorporates gene-level scores from polygenic priority score (PoPS)86, another tool for causal gene prioritization. PoPS prioritizes genes by identifying gene features that are enriched for genetic association across the entire GWAS—such as pathway memberships, protein–protein interaction networks, co-expression data and expression in various tissues and cell types—then scoring each gene based on which of these enriched features it has. PoPS identifies these trait-relevant gene features by leveraging multi-marker analysis of genomic annotation87 to aggregate variant-level GWAS P values into gene-level association scores, then using a regression model to identify which gene features are predictive of high multi-marker analysis of genomic annotation scores.
FLAMES integrates these two evidence streams by multiplying the gradient boosting scores with scaled PoPS scores, effectively upweighting genes supported by both locus-specific and GWAS-wide evidence. At the recommended confidence threshold (scaled FLAMES score >0.248), FLAMES nominated a candidate causal gene at 19 of the 26 loci from the primary meta-analysis, while abstaining from prediction at the remaining 7 loci owing to a lack of confidence in the causal gene.
deCODE pipeline
As a second strategy to nominate candidate causal genes, we performed functional annotation and QTL analyses on all 26 lead variants from the primary meta-analysis, as well as all variants in high LD (r2 ≥ 0.8 and within ±1 Mb) with these lead variants.
We used variant effect predictor88 to attribute to the studied variants the most severe predicted variant consequences in canonical and non-canonical transcripts. We classified as high-impact variants those predicted as start-lost, stop-gained, stop-lost, splice-donor, splice-acceptor or frameshift, collectively called loss-of-function variants, whereas moderate impact variants are those predicted to otherwise affect coding or splicing of a protein (missense).
For all lead and correlated variants, we studied their association with (1) mRNA expression (top local expression quantitative trait locus (eQTL), splicing QTL or alternative polyadenylation QTL) in multiple tissues analyzed at deCODE, in addition to data from GTEx89 and other public datasets, and (2) plasma protein levels (top cis-protein QTL) identified in large proteomic datasets from Iceland and the UK90. RNA sequencing was performed on whole blood from 17,848 Icelanders and on subcutaneous adipose tissue from 769 Icelanders, respectively. Gene expression was computed based on personalized transcript abundances using kallisto91. The association between sequence variants and gene expression (cis-eQTL) was tested via linear regression, assuming additive genetic effect and normal quantile gene expression estimates, adjusting for measurements of sequencing artifacts, demographic variables, blood composition and PCs92. The gene expression PCs were computed per chromosome using a leave-one-chromosome-out method. All variants within 1 Mb of each gene were tested.
The Icelandic proteomics data were analyzed using the SomaLogic SOMAscan v4 proteomics assay that scans 4,907 aptamers, measuring 4,719 proteins in samples from 35,892 Icelanders with genetic information available at deCODE Genetics90,93. Plasma protein levels were standardized and adjusted for year of birth, sex and year of sample collection (2000–2019)90,93. The UK proteomics dataset was analyzed using the Olink Explore 3072 proximity extension assay platform with 2,941 immunoassays characterizing 2,925 proteins in 54,265 participants in the UK Biobank90,93.
Phenome-wide associations of lead variants
For each of our 26 lead variants from the primary meta-analysis, we looked up genome-wide significant (P < 5 × 10−8) associations in the GWAS Catalog that overlapped either the variant itself or any proxy variants in high LD (r2 ≥ 0.8 and within ±1 Mb). To account for the GWAS Catalog containing large numbers of closely related phenotypes, we grouped related traits for visualization. Specifically, we removed GWAS results obtained through MTAG analysis or performed on multiple traits combined, identified by ‘and’ or ‘or’ in their trait names. We then standardized trait names by removing words in parentheses and laterality descriptors (‘left’ and ‘right’), and standardized the spelling of ‘neutrophill’ to ‘neurophil’. Finally, we grouped trait names that became identical after this standardization.
We also looked up the P values of our 26 lead variants in preexisting GWAS summary statistics for each of 330 disease endpoints from the November 2024 release of MVP–Finngen–UKBB (https://mvp-ukbb.finngen.fi/about), a GWAS meta-analysis of Million Veteran Program, FinnGen and the UK Biobank. We performed Bonferroni correction across the 330 diseases and 26 variants tested, leading to a significance threshold of P = 0.05/8,580 ≈ 5.83 × 10−6.
We applied the same approach to 124 GWAS of drug prescription endpoints from FinnGen’s DF12 release, in which each GWAS was on whether an individual was ever prescribed any drug from a particular category (for example, analgesics). We performed Bonferroni correction across the 124 drug categories and 26 variants tested, leading to a significance threshold of P = 0.05/3,224 ≈ 1.55 × 10−5.
Heritability
For the European-only meta-analysis, we estimated single-nucleotide polymorphism-based heritability—the proportion of phenotypic variance explained by the aggregated effect of single-nucleotide polymorphisms across the autosomal genome—using stratified LD score regression81 with the standard 52-annotation baseline model mentioned above. We used LD scores computed from the European-ancestry subsample of 1000 Genomes Phase 3 (precomputed by LDSC authors), and subset to variants present in Europeans in HapMap 3. The observed-scale heritability was calculated from the LD score regression slope. As with all LD score regression-based analyses in this study, we used the effective sample size Neff rather than the raw N to avoid bias due to case–control imbalance.
Tissue and cell-type enrichment
We used LDSC-SEG94 to infer tissue and cell-type enrichments for the European-only meta-analysis. As for the global heritability analysis, we used LD scores computed from the European-ancestry subsample of 1000 Genomes Phase 3 (precomputed by LDSC authors) and the 52-annotation baseline model, subsetting to variants present in Europeans in HapMap 3. LDSC-SEG is a variant of stratified LD score regression that analyzes each tissue or cell type in turn, adding a binary annotation to the 52-annotation baseline model for each one. This annotation flags variants located within or within 100 kb of the 10% of genes most specifically expressed in that tissue or cell type. The heritability enrichment for each tissue or cell type is computed by dividing its per-SNP heritability (derived from the LD score regression slope) by the average per-SNP heritability across all SNPs in the analysis.
To compute tissue enrichments, we used annotations provided by the LDSC-SEG authors that were derived for each of 53 tissues in the GTEx project, in which the top 10% of specifically expressed genes were defined by ranking genes by their differential expression t-statistic for expression in that tissue versus all other tissues. For GTEx brain tissues, the comparison was instead between that tissue and all non-brain tissues.
To compute cell-type enrichments, we generated gene sets from a comprehensive single-cell transcriptomic atlas of the mouse, PanSci, which profiles over 20 million cells from 14 organs and tissues95. From this dataset, we included 119 cell types and 7 broad lineages that were represented by at least 1,000 cells after quality control (filtering to cells with ≤10% mitochondrial reads, ≥100 genes detected and non-zero Malat1 expression). Our analysis was restricted to a universe of 16,404 protein-coding genes with identifiable human orthologs. For each cell type and lineage, we defined its specifically expressed genes as those with a detection rate >10% within that cell type and a fold change >2 compared with all other cells.
Genetic correlation
We performed genetic correlation between our leave-FinnGen-out European-only fibromyalgia meta-analysis and each of 1,284 clinical endpoints from FinnGen Release 13 with https://www.finngen.fi/en/researchers/clinical-endpoints, using LD score regression96. We calculated SNP heritability with LDSC for all 2,691 endpoints and performed genetic correlation only on those endpoints with significant heritability, resulting in a list of 1,284 endpoints. Variants were aligned and filtered to HapMap release 3 SNPs, and the LD reference was obtained from the European population of the 1000 Genomes Project Phase 3 release. We used the 80th percentile of per-variant effective sample sizes (Neff) to munge the summary statistics.
To avoid diluting the results with low-specificity phenotypes, we then excluded endpoints containing any of the (case-insensitive) keywords ‘other’ (for example, ‘Other disorders of ear’), ‘unspecified’ (for example, ‘Disorder of external ear, unspecified’), classified (for example, ‘Other viral diseases, not elsewhere classified’), ‘any’ (for example, ‘Any mental disorder’) or ‘all’ (for example, ‘All kidney diseases’), as well as derived endpoints for which the listed category did not correspond to a chapter of the ICD. We also excluded fibromyalgia itself as an endpoint. These exclusions reduced the number of endpoints tested from 1,284 to 855, resulting in a Bonferroni significance threshold of P = 0.05/855 ≈ 5.85 × 10−5.
For visualization purposes, we grouped endpoints by their ICD chapter. We grouped four categories with few significant endpoints—‘III Diseases of the blood and blood-forming organs and certain disorders involving the immune mechanism’, ‘XVIII Symptoms, signs and abnormal clinical and laboratory findings, not elsewhere classified’, ‘XXI Factors influencing health status and contact with health services’ and ‘XXII Codes for special purposes’—into an ‘Other’ category. Because fibromyalgia has been previously proposed to share etiology with autoimmune disorders, we also created a separate ‘Autoimmune’ category, as these disorders would otherwise be scattered across ICD chapters based on the affected organ system.
We also estimated the genetic correlation between fibromyalgia and the subtypes of asthma and RA using summary statistics available at deCODE. In line with other LDSC-based analyses, we used per-variant effective sample sizes, filtered variants to HapMap release 3 SNPs and used precomputed LD scores for European populations (downloaded from https://data.broadinstitute.org/alkesgroup/LDSCORE/eur_w_ld_chr.tar.bz2).
PRSs
We computed PRSs for the leave-UK Biobank-out European and multi-ancestry meta-analyses, using PRS-CS97. PRS-CS takes GWAS summary statistics and an LD reference panel as input, and applies Bayesian shrinkage to the variants’ effect sizes to infer posterior effect sizes, which are used as the weights of the PRS. These weights can then be scored on any cohort of interest, which should be non-overlapping with the cohorts used to create the PRS to avoid bias.
We used the European-ancestry subsample of 1000 Genomes Phase 3 as the reference panel, removing indels from the summary statistics before the analysis. We ran PRS-CS with the ‘–n_burnin 5000’ and ‘–n_iter 10000’ options to specify 5,000 burn-in iterations and 10,000 total iterations of the Markov Chain Monte Carlo sampling. We then scored the PRS weights on the UK Biobank cohort with the ‘–score’ option from plink version 2.0.0. As PRS-CS requires a single sample size across all variants, the sample size we provided to PRS-CS was the 80th percentile of per-variant effective sample sizes (Neff), as recommended by others for PRS analyses98.
To assess the accuracy of the PRS across ancestries, we scored the PRS weights from the multi-ancestry leave-UK Biobank-out meta-analysis on each of the 3 largest genetic ancestries from the UK Biobank (European, Central and South Asian, and African). We assessed PRS accuracy via the AUC, as well as by computing odds ratios for each quintile of polygenic risk relative to the middle quintile, and prevalence of fibromyalgia within each quintile.
Statistics and reproducibility
No statistical method was used to predetermine sample size. No data were excluded from the analyses. The experiments were not randomized. The investigators were not blinded to allocation during experiments and outcome assessment.
Reporting summary
Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.
