Genome-wide association study on pharmacological outcomes of musculoskeletal pain in UK Biobank

The pharmacological management of musculoskeletal pain starts with NSAIDs, followed by weak or strong opioids until the pain is under control. However, the treatment outcome is usually unsatisfying due to inter-individual differences. To investigate the genetic component of treatment outcome differences, we performed a genome-wide association study (GWAS) in ~23,000 participants with musculoskeletal pain from the UK Biobank. NSAID vs. opioid users were compared as a reflection of the treatment outcome of NSAIDs. We identified one genome-wide significant hit in chromosome 4 (rs549224715, P = 3.88 × 10−8). Suggestive significant (P < 1 × 10−6) loci were functionally annotated to 18 target genes, including four genes linked to neuropathic pain processes or musculoskeletal development. Pathway and network analyses identified immunity-related processes and a (putative) central role of EGFR. However, this study should be viewed as a first step to elucidate the genetic background of musculoskeletal pain treatment.


INTRODUCTION
Chronic musculoskeletal pain is one of the most frequent causes of suffering and disability [1].The nature of musculoskeletal pain can be nociceptive or neuropathic, for which the corresponding pain management differs.The treatment of nociceptive musculoskeletal pain follows the WHO three-step analgesic ladder [2]: the first treatment step is non-opioid analgesics, such as nonsteroidal anti-inflammatory drugs (NSAIDs); the second step is weak opioids for mild to moderate pain, such as tramadol; the third step is strong opioids for moderate to severe pain, such as morphine.
Unfortunately, effective pain management is challenged by inter-individual differences, with unsatisfied pain treatment rates ranging from 34 to 79% [3].The underlying factors of ineffective pain treatment are multifactorial, including demographic characteristics (age, sex, socioeconomic status) [4,5], lifestyle (smoking and alcohol intake) [6], comorbidities (psychological status) [7], and genetic factors.The genetic background of pain treatment outcomes has been investigated using candidate gene approaches.Some drug-metabolizing genes are associated with treatment outcomes for specific drugs, e.g., CYP2D6 and codeine [8].In addition, genes implicated in pain (sensitivity) may contribute to pain treatment outcomes because greater pain sensitivity is associated with increased opioid use [9] and poorer chronic pain treatment outcomes [10].
However, none of these findings predict pain treatment outcomes sufficiently to optimize pain treatment in a clinical setting.Furthermore, these studies are limited by small gene panels and sample size and report contradictory results [11].The most investigated genetic variant is the 118 A to G base pair change in the OPRM1 gene (rs1799971).Genetic variants in OPRM1 are thought to influence the opioid response by altering the µ-opioid receptor binding affinity of exogenous opioid ligands, such as morphine [12].The G allele was associated with higher opioid dosing [13,14] but was shown to be protective against pain in other studies [15,16].Therefore, definitive conclusions on these genetic associations cannot be drawn yet, and a non-hypothesis driven approach in a large population is needed.Except for several recent, successful large-scale GWASs of chronic pain phenotypes [17][18][19], the number of GWASs focusing on pain treatment outcomes is still limited.Moreover, the most frequently used phenotype in GWASes investigating pain treatment is the response to certain drugs for acute pain (e.g., analgesic requirement or pain relief score after surgery [20,21]), but long-term pain treatment outcomes are less investigated.
This study sought to identify genetic variants associated with switching to a higher analgesic ladder in people with musculoskeletal pain from the UK Biobank.A GWAS was performed including subjects treated according to the WHO analgesic ladders, and comparisons were made between NSAID and opioid users as a reflection of pain treatment outcome.

METHOD
We conducted a GWAS comparing NSAID users and opioid users using data from the UK Biobank, and post-GWAS analyses were performed for suggestively significant (P < 1 × 10 −6 ) signals.

Study population
The UK Biobank is a general population cohort with over 0.5 million participants aged 40-69 recruited across the United Kingdom (UK) [22].Recently released primary care (general practitioners, 'GPs') data provides longitudinal structured diagnosis and prescription data, which were used for phenotype definition.At the time of analysis, the interim release of GP data was available, which contained data on approximately 45% of the UK Biobank participants.UK Biobank obtained informed consent from all participants.

Phenotype definition
Figure 1 describes the phenotype definition.To define patients on NSAID and/or opioid treatment, we first extracted all participants with musculoskeletal pain records and participants with pain prescription records from the GP data (see Table S1 for the diagnosis codes and Table S2 for the pain prescriptions (NSAIDs and opioids) codes included in this study).Then, these two datasets were merged by pseudonymized subject ID, GP data provider, and day.As we focus on long-term pain treatment, participants with more than one musculoskeletal pain diagnosis record were included.Patients were eligible for inclusion in the study when they had a pain prescription record occurring on the same day as the diagnosis to ensure that the prescriptions are indeed for musculoskeletal pain treatment.Nociceptive musculoskeletal pain diagnoses were selected in READ code rather than ICD10 codes because the READ code is used throughout the GP data, while the ICD10 code is only available for hospital inpatient records.Therefore, the READ code was used to link diagnosis data with the prescription data in the GP.
The outcome used for the analysis was defined as a dichotomous score (case/control): NSAID users were defined as controls and opioid users as cases.Opioid users were analyzed as a whole because the strong opioid user group is small (n = 365) and assumes mechanistic similarities between weak and strong opioids.Participants who did not meet the following two quality control (QC) steps were removed.First, participants with only one treatment event were removed to safeguard the inclusion of only participants with relatively long-term treatment.Second, a chronological check was applied for the first prescription of each ladder to ensure that the treatment ladder was correctly followed, i.e., opioids followed initial NSAID use.As the GP data is longitudinal, by using this definition, we could distinguish between patients who stay at NSAID treatment and those switching to the next level of the analgesic ladder.The script for defining analgesic ladder switching phenotype can be found in the GitHub repository (https://github.com/lisongmiller/UKB_GWAS_pain_treatment_outcome).

Genotyping and quality control
Genotyping procedures and PCA-based ancestry inference have been described in detail elsewhere [23].Routine QC steps for genetic markers on autosomes included removal of single nucleotide polymorphisms (SNPs) with (1) an imputation quality score less than 0.8, (2) a minor allele frequency (MAF) less than 0.005, (3) a Hardy-Weinberg equilibrium (HWE) test P-value less than 1 × 10 − 6, and (4) a genotyping call rate less than 0.95.
QC steps for genetic markers on the X chromosome's pseudo autosome region (PAR) were the same as autosomes.For non-PAR of the X chromosome, additional QC steps included setting heterozygous haploid genotypes as missing for males, excluding multi-allelic SNPs, and excluding variants with significantly different MAFs between males and females in the NSAID user group (P < 0.05/#SNPs), and variants that violated HWE were examined in the NSAID user group using only females.
Routine QC steps for the samples include the removal of participants with (1) inconsistent self-reported and genetically determined sex, (2) missing individual genetic data with a frequency of more than 0.1, and (3) putative sex-chromosome aneuploidy.Participants were excluded from the analysis if they were considered outliers due to missing heterozygosity, not white British ancestry based on the genotype, and missing covariate data.

Definition of covariates
The following variables from the UK Biobank data set were used for the covariate definition: (1) depression history, which was defined as "YES" if depression records were found in self-reported, inpatient hospital or GP data, and (2) drinking frequency, which was derived from data field 1558: "Daily or almost daily" or "Three or four times a week" was defined as high drinking frequency, other values except for "Prefer not to answer" were defined as low drinking frequency.Differences in categorical covariates were tested by a χ² test.Differences in continuous covariates were compared by t-test in SAS version 9.4 (SAS Institute Inc., Cary, NC).

Genome-wide association study
A GWAS was conducted using binary phenotypes, i.e., NSAID users (controls) versus opioid users (cases), we will refer to this analysis as the primary analysis.For markers on the autosomal chromosomes and PAR region of the X chromosome, GWAS on switching to a higher analgesic ladder was conducted using a linear function in GCTA [24], adjusting for age, sex, BMI, depression history, smoking status, drinking frequency, assessment center, genotyping array, and the first ten principal components (PCs).The selection of PCs was based on a scree plot (Fig. S1).
To examine the nature of pain between groups, all nociceptive musculoskeletal pain diagnosis codes were grouped into one of the following categories: inflammatory, mechanical, and mechanism not specified.The percentage of subjects in each diagnosis category was compared by a χ² test.
Markers on the non-PAR region of the X chromosome were analyzed by a sex-stratified analysis in XWAS [25].A p-value less than 5 × 10 considered genome-wide significant, and P-values between 1 × 10 −6 and 5 × 10 −8 were defined as suggestively significant.FUMA was used to define the lead SNPs (SNPs with the smallest P-value in each locus) and independent significant signals (SNPs in LD (r2 > 0.6) with lead SNPs and remaining significant after conditioning on the lead SNPs in each locus).

Sensitivity analysis and robustness analysis
To test the effect of adjusting study-specific covariates to the association results, a sensitivity analysis was performed by removing all the studyspecific covariates.Therefore, only the following covariates were included in a sensitivity analysis: sex, age, 10PCs, array type, and assessment center.
To validate the GWAS results, a robustness analysis was performed by splitting the sample randomly into two equally sized subsets five times for the lead SNPs (n = 8) as described by Janita Bralten et al. [26].In addition, a more stringent P-value threshold was obtained by a permutation test.In each permutation, genetic associations were performed based on randomly shuffled phenotypes and adjusted for the same covariates as in the primary analysis, and the lowest p-value was recorded.The permutation was run 5000 times, and the permuted threshold is the value at the 5th percentile of the distribution.

Functional annotation
Bayesian fine-mapping of lead loci.As GWAS identified lead variants are not always the causal variants affecting the trait, the true causal variant can be the SNPs correlated to the lead variants through linkage disequilibrium (LD).To identify this, Bayesian fine-mapping in PAINTOR [27] was performed which leverages both the association strength and genomic functional location to prioritize causal variants.Lead SNPs were analyzed to identify the most likely causal SNPs in each locus.PAINTOR calculates the posterior probability (PP) of causality for SNPs in each genomic region by leveraging the strength of association (Z score) and the LD pattern.The 1000 Genomes (Phase 3) were used for LD matrix calculation.The calculated PP for each SNP was sorted from high to low, and variants together reaching a PP of 0.95 were used to define 95% credible sets.
Functional annotation of SNPs present in the 95% credible sets.To evaluate the impact of the non-coding variants in the 95% credible sets, HaploReg v4.1 [28] was used to annotate these SNPs for regulatory functions.Specifically, the analyzed regulatory functions were (1) the presence of exonic, nonsynonymous variants in high LD (r 2 ≥ 0.8), (2) overlap with epigenetic histone marks of active enhancers (H3K4me1 and H3K27ac) and active promoters (H3K4me3 and H3K9ac), and (3) the sensitivity to DNase.As histone marker overlap is tissue-specific, relevant cell lines were selected from the complete data set (see Table S3).Besides regulatory functions, potential pleiotropy effects (previously reported associations with other phenotypes) of the variants were investigated in Haploreg.For SNPs not available in Haploreg, proxy SNPs were obtained by LD proxy (https://ldlink.nci.nih.gov/).For loci containing more than ten likely causal variants, only the lead SNP and SNPs with the maximum posterior probability (PPmax) of the SNPs in one locus were annotated.
Gene mapping.To annotate suggestively significant GWAS SNPs (Pvalue < 1 × 10 −6 ) and SNPs in LD (LD > 0.6) with them to genes, a bioinformatic tool FUMA was used.Three strategies were adopted in FUMA: positional mapping, expression quantitative trait loci (eQTL) mapping, and chromatin interaction mapping.For the positional mapping, SNPs were functionally annotated to known protein-coding genes based on physical distance (within a 10 kb window).For eQTL mapping, SNPs were functionally annotated to genes up to 1 Mb away based on known cis-eQTLs.As gene expression is tissue-specific, we selected the following tissues for mapping: brain, muscle, kidney, liver, nerve, skin, and fibroblast.Significant eQTLs were defined as eQTLs with a false discovery rate (FDR) < 0.05.Finally, chromatin interactions for gene mapping were assessed.Chromatin interaction can occur in two genomic regions that are spatially close when DNA folds together, even if the genomic regions are at a long-range physical distance.Genes in regions of chromatin interaction containing candidate SNPs were assessed in the same tissues as the eQTL mapping.An FDR < 1 × 10 −6 was defined as a significant interaction.It is noteworthy that FUMA links SNPs to genes by combining information from biological data sources to generate hypotheses for further functional validation experiments.Therefore, annotated genes are not guaranteed to be causally linked to the investigated SNP.
Heritability and power calculation.Narrow-sense heritability was calculated by the GREML method in GCTA, which first calculates the genetic relationship matrix from all the autosomal SNPs, then performs a REML (restricted maximum likelihood) analysis [29,30].The same covariates used in the GWAS analysis were included in the heritability analysis.The power of heritability was calculated using GCTA-GREML Power Calculator (https:// shiny.cnsgenomics.com/gctaPower/)assuming a disease risk in the population of 0.5, trait heritability as 0.08, and α level as 0.05.
The power of GWAS was calculated using the CaTS power calculator (http://csg.sph.umich.edu/abecasis/cats),assuming an additive model with the following input parameters: a significance level of 5 × 10 −8 , the prevalence of phenotype (opioid use) of 50%, and a relative genotypic risk of 1.135, based on 11,089 cases and 12,726 controls.
Pathway enrichment analysis.To investigate if the genes identified in the GWAS could be linked to biological pathways and networks involved in pain (treatment), a gene-based functional pathways enrichment was performed using Ingenuity Pathway Analysis software (IPA®, QIAGEN Inc., Redwood City, California, U.S.).IPA is based on prior knowledge of direct and indirect gene relationships from experimentally observed data in mammals and all cell types.A gene-based P-value was computed twice using the gene analysis function in MAGMA v1.08.Firstly, a gene-based P-value was calculated using nominally significant SNPs (P-value < 0.05) in the protein-coding region of genes without flanking regions.Secondly, the same analysis was performed, including only nominally significant SNPs in the protein-coding gene regions and 100 kilobases (kb) pair upstream and downstream flanking regions, to take cis-eQTL effects into account [31][32][33].Both analyses were combined, and the smallest P-value of each gene was selected for pathway and network analyses with IPA (see gene list in Table S4).Pathways with an FDR < 0.05 were considered statistically significantly enriched.To illustrate the core molecules in the networks, a radial plot was generated of the merged top five networks.
Genetic correlation analysis.Genetic correlation is the proportion of variance that two traits share due to genetic causes, i.e., the correlation between the genetic influences on a trait and the genetic influences on a different trait.Identifying genetic correlations can estimate the degree of pleiotropy or causal overlap between complex traits and diseases [34].
Genetic correlations between switching to a higher analgesic ladder and other complex traits were investigated by linkage disequilibrium score regression through LD Hub v1.9.3.The tested traits were selected from the LD hub, and the following categories were selected: education, anthropometric, sleeping, psychiatric, personality, cognitive, autoimmune, smoking behavior, kidney, neurological, and UKB phenotypes.Correlations with P-values less than 8.4 × 10 −5 (0.05/596) were considered significant.Since the top nominally significant correlations were overrepresented by pain (category 100048, and data field 6154 in the UKB) and socioeconomic status phenotypes, the percentage of nominally significant (P < 0.05) correlations in these two categories were compared with all other categories by a χ² test.The socioeconomic status phenotypes consisted of qualifications (data-field 6138), employment (category 100064) in the UKB, and all education phenotypes in the LD hub [35][36][37][38].
Ordinal GWAS.Besides the binary case/control analysis, an additional GWAS was performed using an ordinal outcome.For the 'ordinal phenotype', an ordinal score of '1', '2', or '3' was assigned to NSAID users (persons only using NSAIDs), weak opioid users (persons using NSAIDs and weak opioids), and strong opioid users (persons using NSAID, weak opioid, and strong opioids), respectively.Also for the ordinal analysis, patients with one treatment event and not following chronological treatment were removed.The ordinal regression analysis was conducted in OrdinalGWAS [39] using the same covariates as the binary analysis.
Subtype GWAS and secondary GWAS.As inflammatory pain is an important subtype of pain, a subtype GWAS was carried out for this phenotype specifically.Only inflammatory nociceptive musculoskeletal pain diagnosis codes were used for participant selection.Some participants with inflammatory nociceptive musculoskeletal pain received pain treatment for other types of pain (e.g., mechanical nociceptive musculoskeletal pain or not specified) over time.These participants were excluded from the analysis.
Moreover, as diagnostic codes were often not repeatedly recorded [40] or reported as repeat prescriptions, a secondary GWAS was performed for pain medication users with less strict diagnosis criteria than the primary analysis.To ensure a relatively homogenous population, we included participants with at least one nociceptive musculoskeletal pain treatment event but without any other limitations on their pain treatment purposes.To focus on long-term treatment and remove outliers, participants were removed if they only had one prescription record or prescription record numbers on a log scale outside the 1.5 inter-quantile range.In addition, participants already included in the primary GWAS analysis were not included in this analysis.
Both the subtype and secondary GWASs were performed using the binary phenotype and following the same procedure as the primary analysis, but only autosomal markers were examined.To investigate whether the identified loci in these two analyses overlapped with the primary analysis, lead SNPs with a suggestive threshold (P < 1 × 10 −6 ) were compared with the primary GWAS signals for LD correlation (r2 > 0.6) in LDpair (https://ldlink.nci.nih.gov).

GWAS
After quality control, we identified 12,726 NSAID users (control) and 11 089 opioid users (cases) in the UK Biobank dataset.Table 1 summarizes the demographics of the cases and controls, and all tested covariates were found to be significantly different (P < 0.0001).
There were 9,435,994 SNPs available for GWAS analysis after quality control.The genomic control value (lambda) was 1.008.One intergenic locus located at chromosome 4 reached genome-wide significance, in which the most significant SNP was rs549224715 (P = 3.92 × 10 −8 ) (Fig. 2, Table 2).Seven loci surpassed the suggestive P-value threshold (P < 1 × 10 −6 ), and no other independent SNPs (SNPs remaining significant after conditioning on lead SNPs in the locus) were identified in each locus.In addition, the ordinal GWAS results were consistent with the GWAS using binary outcomes (Fig. S2).However, we did not find any genome-wide significant loci or overlapped suggestively significant loci with the primary GWAS (Fig. S3, Fig. S4, Table S5, Table S6).
Our study had 80% power to identify SNPs with a risk allele frequency of 5% and a genotypic relative risk of 1.135.The SNP heritability was 0.16 (P-value = 0.16), and the power for this analysis was 96%.

Sensitivity analysis and robustness analysis
In the sensitivity analysis, the lead SNP rs549224715 in the primary GWAS remained the strongest signal, and two more lead SNPs reached genome-wide significance (rs143781228 and rs34147893 (in complete LD with SNP rs13133042 identified in the primary analysis)).Three other lead SNPs in the primary GWAS (rs12694371, rs73062440, and rs10264795) remained suggestively significant, except for the SNPs rs73126966 and rs13218801 (Table 2).
In the robustness analysis, six of eight suggestively loci passed at least the nominal significance threshold (P < 0.05) for all ten iterations.Two lead SNPs (rs13218801 and rs73062440) failed one out of ten iterations (Table S7).
In the permutation analysis, the permutated P-value at the 5th percentile is 1.61 × 10 −8 (Fig. S5).None of the lead SNPs passed this threshold.

Functional annotation
Statistical fine-mapping of loci and functional annotation of SNPs.In five out of eight loci, the lead SNPs in the locus had the maximum PP (PP max ) (Fig. S6).The regulatory function annotation results suggested that most genetic variants were potentially involved in transcriptional regulatory modulation (Fig. S6).
We assessed whether the SNPs in the 95% credible sets were previously reported to be associated with other traits.However, no pleiotropic effects were identified.Gene mapping.A total of 18 unique annotated genes were identified by SNPs in the 95% credible sets (Table 3).Five genes were annotated by genomic location, seven genes were identified by cis-eQTL mapping, ten genes were annotated as SNPs in regions where 3D chromatin interactions occurred, and four genes were identified by at least two mapping strategies.

Pathway enrichment
Pathway enrichment analysis in IPA prioritized 15 significant pathways with an FDR < 0.05, in which the top prioritized pathways were mainly implicated in the immunological response.(Table S8).
The network analysis yielded a total of 25 prioritized networks.The top network contained 33 genes with the EGFR protein in the center.EGFR remained in the center after merging the five networks with the lowest P-value (Table S9, Fig. S7).

DISCUSSION
In this study, we identified one genome-wide significant hit and seven loci with suggestive significance associated with analgesic ladder switching from NSAID to opioid.Although pain or pain treatment is characterized by sex differences, i.e., females are more vulnerable to pain and opioid use [5], no significant signals were found on the X chromosome.
The functional link between the genome-wide significant SNP (rs549224715) on chromosome 4 and switching to a higher pain treatment ladder remains unclear.The nearest gene of this locus, CWH43, is associated with Seckel Syndrome, characterized by growth delays before birth.Another gene, TXK, was annotated by eQTL to this SNP which has previously been linked to rheumatoid arthritis [41].Considering that our investigated phenotype is switching from NSAID to opioid use, this association is worth further validation and investigation.
Two genes identified from other loci are prioritized as they are involved in neuropathic pain.NPTX2 was identified by both eQTL mapping and gene-based analysis with the lowest P-value (2.71 × 10 −5 ) (Table S2).NPTX2 is down regulated in the brain in induced chronic neuropathic pain [42] and induced endometriosis [43] mouse models.The other gene is IPCEF1, is previously related to neuropathic pain conditions [44].In addition, we identified two genes linked to muscular or skeletal dystrophy: SGCB and FN1.These genes are of interest as musculoskeletal dystrophy is characterized by pain.
None of the identified genes were implicated in the metabolism and working mechanisms of specific NSAIDs, which is perhaps not unexpected.As subcategories in NSAIDs, such as non-selective NSAIDs and selective COX-2 inhibitors, were analyzed as a whole, this may have diluted signals related to a specific NSAID.However, stratified analyses per drug were not practical as the groups would become too small to obtain sufficient power.
The narrow-sense heritability result in our study is in agreement with similar phenotypes, opioid response (60% in cold pressorinduced pain and 12% in heat pressor) in a twin study [45] and chronic pain (0.08 to 0.31) [7].Since we have sufficient power for the heritability analysis, the insignificant narrow-sense heritability result suggests that genetic factors might not play an important role in switching to a higher analgesic ladder compared with other well-known environmental, psychological, and socioeconomic factors influencing pain and pain treatment [46].However, our results should be further validated in a larger sample size considering the fact that a large number of common variants with very small effect sizes contribute to the complex traits.Take height as an example, a close heritability estimation of height to the pedigree-based heritability estimation was made possible with larger sample sizes with n > 100,000 for height in a recent GWAS meta-analysis [47].
Our study might add to the current evidence of the biological mechanisms of pain treatment.Most variants in the 95%  Demographics of NSAID users (control) and opioid users (cases) in the UK Biobank.Categorical covariates are represented as frequency (percentage) and compared by the χ2 test.Continuous covariates are presented as mean (standard deviation) and compared by an independent t-test.Depression history was defined as "YES" if depression records were found in the self-reported, inpatient hospital, or GP data.Drinking frequency was defined from data field 1558, "Daily or almost daily" or "Three or four times a week" defined as high drinking frequency, other values except for "Prefer not to answer" defined as low drinking frequency.
a Percentage of subjects within a certain category of pain in each group.Footnote a: one subject could have more than one type of diagnosis, so the percentage sum is not equal to 1.
credible sets showed potential transcription regulatory functions, which aligns with research indicating that epigenetic changes are involved in chronic pain [48] and pain treatment [49].For instance, epigenetic restructuring in a candidate gene (OPRM1) and global DNA methylation was observed after opioid use [50,51].Furthermore, the network analysis identified EGFR, and preliminary studies suggest a role of EGFR in pain modulation in preclinical studies [52] and analgesic effects in clinical settings [53].However, the results from this analysis should be interpreted carefully as the input consisted of nominally significant genes from the GWAS analysis.In addition, the overrepresentation of pain phenotypes in genetic correlation analysis also matched reports that opioid users tend to have more chronic and severe pain conditions [54].
Although six lead SNPs passed the robustness analysis, only the lead SNP passed the permutation threshold in the sensitivity analysis, and none of the SNPs survived permutation analysis which might suggest these results can still be false positive findings.In addition, we failed to replicate our findings in the secondary analysis.Unfortunately, replicating the results in other independent cohorts is difficult due to the limited number of publicly available large-scale data similar to UK Biobank and the lack of cohorts measuring similar outcomes.It might still be worth validating our results in a large cohort with a clear outcomes definition, such as the ongoing Pain Predict Genetics cohort in our center (NCT02383342).
The merit of our study is that we have a large sample size investigating analgesic ladder switching from NSAIDs to opioids.However, the limitation is that we utilized a derived phenotype, so we cannot distinguish whether switching ladders is because of pain progress, poor treatment response to certain drugs, or a combination of both factors.However, it does not matter whether SNPs mapped to candidate genes using FUMA.SNPs in LD (r2 > 0.6) with lead SNPs were mapped to genes.posMapSNPs, the number of SNPs mapped to this gene based on positional mapping; eqtlMapSNPs, the number of SNPs mapped to this gene based on eQTL mapping; eqtlMapminP, the minimum eQTL P-value of mapped SNPs; eqtlMapminQ, the minimum eQTL FDR of mapped SNPs; ciMap, "Yes" if the gene is mapped by chromatin interaction mapping; minGwasP, the minimum P-value of mapped SNPs.Overview of the lead SNPs passing suggestive significance in the primary GWAS for NSAID versus opioid users.Bold font indicates the SNP that passed the genome-wide significant threshold (5 × 10 − 8).CHR:POS physical position of the SNP, A1 effect allele, AF1 effect allele frequency, BETA (SE) effect size of SNP, and standard error (SE).# Identified SNP in the sensitivity analysis is rs34147893, which is in complete linkage disequilibrium with rs13133042.
the genes might reflect pain severity or pain treatment outcome, as they still have the potential to predict analgesic ladder switching.Despite the limitations in phenotype definition, the group characteristics are similar to previous publications, with a roughly even share of NSAID users and opioid users in the population [55], and the reported risk factors for using opioids are also in line with previous literature [54,55].
In conclusion, we identified one locus achieving genome-wide significance for a derived pain treatment outcome phenotype.Some identified genes could be linked to neuropathic pain and musculoskeletal development.This study should be viewed as an initial stepping stone for future research.We show a small genetic contribution to analgesic ladder switching.For future research, it might be better to focus on a specific disease or outcome for a specific treatment.

Fig. 2 Q
Fig. 2 Q-Q plot and Manhattan plot of primary analysis for compating NSAID versus opioid users.a Q-Q plot of the GWAS results.The red line indicates the distribution under the null hypothesis, and the shaded area indicates the 95% confidence band.b Manhattan plot of the GWAS results.The red line corresponds to the genome-wide significance threshold of 5 × 10 −8 , whereas the blue indicates the suggestive threshold of 1 × 10 −6 .Lead variants are highlighted as orange diamonds.Variants in one locus (within 400 Kb) are highlighted in orange.

Table 1 .
Demographic of NSAID users (control) and opioid users (cases) in the UK biobank.

Table 2 .
Lead SNPs passing the suggestive significance level in the primary GWAS comparing NSAID and opioid users.