Protocol for finding genetic variation associated with unmeasured traits through GenomicSEM common-factor GWAS.

REAGENT or RESOURCE SOURCE IDENTIFIER Deposited data HapMap3 SNP set The International HapMap3 Consortium 5 https://utexas.app.box.com/s/vkd36n197m8klbaio3yzoxsee6sxo11v GWAS Catalog Cerezo et al. 6 https://www.ebi.ac.uk/gwas/downloads/summary-statistics Software and algorithms ldsc Bulik-Sullivan et al. 3 , 4 https://github.com/bulik/ldsc GenomicSEM Grotzinger et al. 2 https://github.com/GenomicSEM/GenomicSEM MungeSumstats Murphy et al. 7 https://github.com/neurogenomics/MungeSumstats Other Linux operating system NA NA Yale McCleary HPC NA https://docs.ycrc.yale.edu/clusters/mccleary/ This protocol describes how to carry out a common-factor GWAS analysis using GenomicSEM as performed in the associated publication. 1

K.J.A.J. prepared scripts and first draft of the manuscript; all authors reviewed and edited the manuscript.

GenomicSEM is a multivariate method for analyzing complex-trait genetic architecture with a wide range of different applications. 2 Here we describe the common-factor GWAS, one specific application of GenomicSEM, as being applied to find SNPs significantly associated with nociplastic pain, an unmeasured common latent factor that we hypothesize to underlie COPCs. GenomicSEM common-factor GWAS has been generally described and there is also an extensive GitHub wiki page available. GenomicSEM in general is comprised of two steps: first, multivariable LD score regression 3 , 4 (LDSR) is used to estimate a genetic covariance matrix. Off-diagonal elements of this sampling matrix indicate the degree to which sampling errors are associated (i.e., degree to which traits overlap). Then, a structural equation model (path diagram/ model) is specified. Parameters of this model are then estimated through minimizing the difference between the empirical genetic covariance matrix and the matrix implied by the structural equation model. This is a flexible framework with a wide range of applications, and in a common-factor GWAS analysis the two-step procedure is extended, expanding the genetic covariance and sampling matrices to include individual SNP effects. Commonly used abbreviations: GWAS = genome wide association study, COPC = chronic overlapping pain condition, LDSR = linkage disequilibrium score regression, SNP = single nucleotide polymorphism, HPC = high performance computing, HapMap = International HapMap Project haplotype map datasets.

The expected outcomes of this analysis are a set of summary statistics that look very similar to standard GWAS output, and include columns ‘Est’ (estimate, i.e., analogous to a GWAS regression beta value column), and its associated p value. Additional columns in the common-factor GWAS output file include Q (heterogeneity), a measure of whether variants relationship with the underlying latent factor is truly multivariate or is likely operating largely through one specific contributing GWAS trait, p value associated with Q, a column indicating whether the model failed to converge for that SNP (fail), and a note of any warnings generated during that model fitting (‘warning’ column) ( Figure 2 ). Figure 2 Example of GenomicSEM ‘commonfactor’ output Example of GenomicSEM ‘commonfactor’ output

Further information and requests for resources should be directed to and will be fulfilled by the lead contact, Keira Johnston ( [email protected] ). Technical questions on executing this protocol should be directed to and will be answered by the technical contact, Keira Johnston ( [email protected] ). No newly generated materials as a result of this protocol. Packages required for this protocol are freely available (see key resources table ), data generated as part of the primary study (common-factor GWAS and TWAS output) is available, and the accession number for this data is Zenodo: 8117583 . Code examples for each step of our common-factor GWAS with GenomicSEM are provided in a GitHub repository: https://github.com/kjaj/gsem_commonFactorGWAS_example .

GenomicSEM is a flexible framework with a wide range of applications beyond common-factor GWAS. Even within common-factor GWAS, it is usually time-intensive and potentially impossible to assess every potential model manually. While GenomicSEM can be applied in a hypothesis-free, exploratory fashion, it is still important to have background knowledge of the traits included in analyses, which can guide the decisions in deciding on the common-factor GWAS model to proceed with. For example, one giant common-factor model of every complex trait for which there is available GWAS is not likely to be useful or informative, and assessing model fit for every possible combination of this large range of traits is not usually feasible. There is significant risk of overfitting – in fitting a large number of models, the best-fit model may not actually generalize well to new data. Solutions to this include comparison of a small number of pre-decided different SEMs, or a-priori choice of one particular model, amongst other approaches. 13 , 14 Another major limitation can be GWAS summary statistics – the authors of the GenomicSEM package discuss in detail sourcing, processing, and understanding GWAS datasets before using them in GenomicSEM analyses, and serious errors can result when this is not carefully considered.

Timing: variable Timing: variable Here we describe setting up the computing environment and software for GenomicSEM analyses. 1. Install GenomicSEM and all dependencies (see key resources table ). 2. Download reference datasets ( key resources table ). Note: For initial sumstats processing and munging, and multivariable LDSR, it is possible to use an interactive job environment on the terminal with R and/or bash. For the ‘sumstats’ function implementation running the model without SNPs, and running the common factor GWAS, this may take sufficient walltime (scale of several hours or more) to require sending jobs to an HPC queue. Depending on analysis size and computing power, it may be possible to carry out all steps on an individual computer, or in an interactive HPC environment running R without scheduling or needing additional resources. For example, we use Yale’s McCleary HPC cluster, a Linux operating system using Slurm for scheduling. For interactive jobs 20G memory was requested and was adequate, for scheduled jobs 28 CPUs 1200G ensured per-chromosome common-factor GWAS with individual SNP effects was limited to ∼1 h walltime. Install GenomicSEM and all dependencies (see key resources table ). Download reference datasets ( key resources table ). Note: For initial sumstats processing and munging, and multivariable LDSR, it is possible to use an interactive job environment on the terminal with R and/or bash. For the ‘sumstats’ function implementation running the model without SNPs, and running the common factor GWAS, this may take sufficient walltime (scale of several hours or more) to require sending jobs to an HPC queue. Depending on analysis size and computing power, it may be possible to carry out all steps on an individual computer, or in an interactive HPC environment running R without scheduling or needing additional resources. For example, we use Yale’s McCleary HPC cluster, a Linux operating system using Slurm for scheduling. For interactive jobs 20G memory was requested and was adequate, for scheduled jobs 28 CPUs 1200G ensured per-chromosome common-factor GWAS with individual SNP effects was limited to ∼1 h walltime. Timing: variable (varies depending on availability of sumstats, necessity of application to consortia for sumstats, and how responsive original GWAS authors may be) Timing: variable (varies depending on availability of sumstats, necessity of application to consortia for sumstats, and how responsive original GWAS authors may be) In this step we describe how to identify GWAS for analyses and how to obtain summary statistics through publicly available download or contact with original study authors. 3. Identify traits of interest and find well-powered GWAS. Considerations include sample size and genetic ancestry of study sample, and type of genotyping chip (see troubleshooting ). 4. Download summary statistics from publicly available sources, or request them from consortia or directly from GWAS study authors. Note: Depending on availability of summary statistics and application requirements or author responsiveness, this step may take months. Identify traits of interest and find well-powered GWAS. Considerations include sample size and genetic ancestry of study sample, and type of genotyping chip (see troubleshooting ). Download summary statistics from publicly available sources, or request them from consortia or directly from GWAS study authors. Note: Depending on availability of summary statistics and application requirements or author responsiveness, this step may take months. Timing: approximately 1 h max Timing: approximately 1 h max In this step we describe how to perform initial quality control checks on summary statistic files. 5. Perform basic checks ensuring files contain a sensible number of SNPs, that all expected columns are present, and that checksums if provided are correct. 6. Assess whether the summary statistics contain rsIDs as SNP variant identifiers. If not, use available maps and/or tools such as the R package ‘MungeSumstats’ to convert SNP IDs to rsID. Note: It is important to know the genome build used in the original GWAS analysis. Sometimes this will be apparent from the sumstats file name, a column in the sumstats, the original GWAS paper, or can be found by spot-checking an rsID position and cross-referencing with dbSNP. If it is still not clear which genome build has been used, this may require reaching out to the original GWAS authors for clarification, and/or close reading of paper(s) associated with the original biobank sample (e.g., if the GWAS was carried out using UK Biobank, information on genotyping can be found in Bycroft et al. 8 ). Finally, there are also R packages (e.g., MungeSumstats) with functions available to compare rsID:CHR:BP in sumstats to HapMap3 GRCh37 and GRCh38 positions to determine sumstat genome build, or this can be carried out manually. 7. Determine if effect size (beta, odds ratio), standard error of effect size, and p value columns are available and that they are they labeled in a ‘standard’ way (e.g., BETA, Beta, beta, or effect). 8. Ensure a minor allele frequency (MAF) column is present. Note: In some GWAS summary statistics only an effect allele frequency (EAF) value is given (i.e., the frequency 0-1 of the allele for which association with the trait is measured in GWAS). This is not necessarily the minor allele (minor allele is defined as the second most common allele at a SNP). One way to assess this is checking if any value in this column is > 0.5, if so, this value is not referring to the minor allele frequency (the minor allele of a biallelic SNP by definition has to be present less than half of the time). MAF can be calculated as 1 – EAF, where EAF > 0.5. 9. Find sample sizes for GWAS and calculate effective sample size if required. Note: This may involve careful reading of the original GWAS publication including supplemental information. If the GWAS summary statistics come from a GWAS meta-analysis of a case-control phenotype (i.e., multiple separate GWAS were performed for the trait IBS (case-control), and then the final sumstats provided and used in your GenomicSEM analyses are a meta-analysis of those multiple separate GWAS), an ‘effective sample size’ ( N e f f ) should be calculated (see github wiki article on effective sample size 9 ). This is necessary because an effective sample size derived from total cases and total controls across each sample population included in the GWAS meta-analysis will be biased by varying ascertainment across samples. 10 An unbiased (effective) sample size estimate is needed in order to calculate an unbiased liability scale SNP-heritability value, which is a key component of all LD-score regression, and a key part of producing the genetic covariance matrix needed to estimate model parameters and perform common-factor GWAS (or any analyses) with GenomicSEM. For example for IBS, we calculated a per-SNP N e f f as: (Equation 1) N e f f = 4 2 ∗ M A F ∗ ( 1 − M A F ) ∗ S E 2 Perform basic checks ensuring files contain a sensible number of SNPs, that all expected columns are present, and that checksums if provided are correct. Assess whether the summary statistics contain rsIDs as SNP variant identifiers. If not, use available maps and/or tools such as the R package ‘MungeSumstats’ to convert SNP IDs to rsID. Note: It is important to know the genome build used in the original GWAS analysis. Sometimes this will be apparent from the sumstats file name, a column in the sumstats, the original GWAS paper, or can be found by spot-checking an rsID position and cross-referencing with dbSNP. If it is still not clear which genome build has been used, this may require reaching out to the original GWAS authors for clarification, and/or close reading of paper(s) associated with the original biobank sample (e.g., if the GWAS was carried out using UK Biobank, information on genotyping can be found in Bycroft et al. 8 ). Finally, there are also R packages (e.g., MungeSumstats) with functions available to compare rsID:CHR:BP in sumstats to HapMap3 GRCh37 and GRCh38 positions to determine sumstat genome build, or this can be carried out manually. Determine if effect size (beta, odds ratio), standard error of effect size, and p value columns are available and that they are they labeled in a ‘standard’ way (e.g., BETA, Beta, beta, or effect). Ensure a minor allele frequency (MAF) column is present. Note: In some GWAS summary statistics only an effect allele frequency (EAF) value is given (i.e., the frequency 0-1 of the allele for which association with the trait is measured in GWAS). This is not necessarily the minor allele (minor allele is defined as the second most common allele at a SNP). One way to assess this is checking if any value in this column is > 0.5, if so, this value is not referring to the minor allele frequency (the minor allele of a biallelic SNP by definition has to be present less than half of the time). MAF can be calculated as 1 – EAF, where EAF > 0.5. Find sample sizes for GWAS and calculate effective sample size if required. Note: This may involve careful reading of the original GWAS publication including supplemental information. If the GWAS summary statistics come from a GWAS meta-analysis of a case-control phenotype (i.e., multiple separate GWAS were performed for the trait IBS (case-control), and then the final sumstats provided and used in your GenomicSEM analyses are a meta-analysis of those multiple separate GWAS), an ‘effective sample size’ ( N e f f ) should be calculated (see github wiki article on effective sample size 9 ). This is necessary because an effective sample size derived from total cases and total controls across each sample population included in the GWAS meta-analysis will be biased by varying ascertainment across samples. 10 An unbiased (effective) sample size estimate is needed in order to calculate an unbiased liability scale SNP-heritability value, which is a key component of all LD-score regression, and a key part of producing the genetic covariance matrix needed to estimate model parameters and perform common-factor GWAS (or any analyses) with GenomicSEM. For example for IBS, we calculated a per-SNP N e f f as: We then calculated total N e f f from the overall Ncase, Ncontrol as: (Equation 2) N e f f t o t a l = 4 ∗ v ∗ ( 1 − v ) ∗ ( N c a s e + N c o n t r o l ) Where v = case proportion. Using N e f f t o t a l to cap each per-SNP N e f f value (max = 1.1 x N e f f t o t a l , min = 0.5 x N e f f t o t a l ). These are capped in this way following recommendations from the GenomicSEM wiki, 9 which are based on Prive et al 11 and minimize the effects of outliers on the final N e f f estimate. To obtain a final single N e f f value per trait (for use in munging step) we take the mean of per-SNP N e f f . For low back pain, we calculated per-cohort N e f f using Equation 2 and took the final N e f f value for the trait to be the sum of these per-cohort effective sample sizes. Further options to obtain N e f f are discussed in the GenomicSEM wiki. For example: require(data.table) dat<-fread("chronic_pain_condition_sumstats_misc/IBS_GCST90016564_buildGRCh37.tsv", data.table=F) dat$MAF 0.5)] 0.5)]) dat$Neff<-4/((2∗dat$MAF∗(1-dat$MAF))∗dat$standard_errorˆ2) Ncases<-53400 Ncontrols<-433201 v<-Ncases/(Ncases+Ncontrols) TotalNeff<-4∗v∗(1-v)∗(Ncases+Ncontrols) #cap at 1.1 of total effective N dat$Neff 1.1∗TotalNeff, 1.1∗TotalNeff, dat$Neff) #lower limit of 0.5 of total effective N dat$Neff<-ifelse(dat$Neff < 0.5∗TotalNeff, 0.5∗TotalNeff, dat$Neff) dat$effect_allele_frequency<-NULL dat$N_CASE<-NULL dat$N_control<-NULL 10. Find sample prevalence for each GWAS trait – for non-meta-analysis GWAS (i.e., those where you did not have to calculate sum of effective sample size), this is simply Ncase/Ntotal. For meta-analysis GWAS you will set this to be 0.5. For continuous traits, set the value to NA. 11. Find population prevalence for each GWAS trait. Note: This value can be obtained per trait from review of the literature, ideally finding large epidemiological studies that assess populations matching GWAS cohorts (in terms of demographics, genetic ancestry) as much as is possible, while maximizing sample size. CRITICAL: Prevalence is distinct from incidence, the former being the proportion of people with the disease or trait in a population at any one time, and incidence being a measure of new cases emerging during a given time period. 12. Save the sumstats files. Note: First, save each set of sumstats in a suitable format (tab-separated tables), one file per trait – these will be input for the multivariable LDSR step. Separately, save another set of tables (1 per chromosome per GWAS trait) – these will be the input for the ‘sumstats’ function step in GenomicSEM. Preparing/ editing GWAS summary statistic files can be carried out on the command line or using R, or a mixture of both. require(data.table) dat<-fread("chronic_pain_condition_sumstats_misc/IBS_GCST90016564_buildGRCh37.tsv", data.table=F) dat$MAF 0.5)] 0.5)]) dat$Neff<-4/((2∗dat$MAF∗(1-dat$MAF))∗dat$standard_errorˆ2) Ncases<-53400 Ncontrols<-433201 v<-Ncases/(Ncases+Ncontrols) TotalNeff<-4∗v∗(1-v)∗(Ncases+Ncontrols) #cap at 1.1 of total effective N dat$Neff 1.1∗TotalNeff, 1.1∗TotalNeff, dat$Neff) #lower limit of 0.5 of total effective N dat$Neff<-ifelse(dat$Neff < 0.5∗TotalNeff, 0.5∗TotalNeff, dat$Neff) dat$effect_allele_frequency<-NULL dat$N_CASE<-NULL dat$N_control<-NULL Find sample prevalence for each GWAS trait – for non-meta-analysis GWAS (i.e., those where you did not have to calculate sum of effective sample size), this is simply Ncase/Ntotal. For meta-analysis GWAS you will set this to be 0.5. For continuous traits, set the value to NA. Find population prevalence for each GWAS trait. Note: This value can be obtained per trait from review of the literature, ideally finding large epidemiological studies that assess populations matching GWAS cohorts (in terms of demographics, genetic ancestry) as much as is possible, while maximizing sample size. CRITICAL: Prevalence is distinct from incidence, the former being the proportion of people with the disease or trait in a population at any one time, and incidence being a measure of new cases emerging during a given time period. Save the sumstats files. Note: First, save each set of sumstats in a suitable format (tab-separated tables), one file per trait – these will be input for the multivariable LDSR step. Separately, save another set of tables (1 per chromosome per GWAS trait) – these will be the input for the ‘sumstats’ function step in GenomicSEM. Preparing/ editing GWAS summary statistic files can be carried out on the command line or using R, or a mixture of both. Timing: 1 h Timing: 1 h In this step we describe how to munge summary statistics and perform multivariable LD score regression. ‘Munging’ is a computing and data science term, is commonly used by ldsc and GenomicSEM authors, and here means ldsc-standard data manipulation and filtering. 13. Take the summary statistic files prepared in the previous step (1 file per GWAS trait, i.e., NOT the sumstats files that are split by chromosome) and use the GenomicSEM R package ‘munge’ function to munge sumstats. Note: Take care to provide the Neff (see 9 ) to the vector of sample sizes, and to maintain trait order (see troubleshooting ). This step can also be completed using ‘munge’ as part of command-line ldsc. The output will be one file per trait, naming convention ‘trait.sumstats.gz’, a compressed tab-separated file containing 1 row per each SNP that meets MAF > 0.01 and INFO (imputation quality score, if available in the summary statistic file input) > 0.9 and columns SNP (rsID) BETA A1 (effect allele) A2 and P. This will be the input for multivariable LDSR ( Figure 1 ). 14. Taking the files produced in step 4 (ending .gz), use the GenomicSEM R package to carry out multivariable LDSR on all the traits you intend to include in the common-factor GWAS. Multivariable LDSR is a means to calculate genetic correlation between all traits, and produces the empirical genetic covariance matrix used for estimating parameters of genomic SEMs. Note: Be sure to assign the sample and population prevalence values (10 & 11), and to maintain trait order and order of corresponding attributes (sample prevalence, population prevalence, trait names). Figure 1 Example of files for GenomicSEM ‘sumstats’ function and files for input in GenomicSEM or ldsc ‘munge’ Take the summary statistic files prepared in the previous step (1 file per GWAS trait, i.e., NOT the sumstats files that are split by chromosome) and use the GenomicSEM R package ‘munge’ function to munge sumstats. Note: Take care to provide the Neff (see 9 ) to the vector of sample sizes, and to maintain trait order (see troubleshooting ). This step can also be completed using ‘munge’ as part of command-line ldsc. The output will be one file per trait, naming convention ‘trait.sumstats.gz’, a compressed tab-separated file containing 1 row per each SNP that meets MAF > 0.01 and INFO (imputation quality score, if available in the summary statistic file input) > 0.9 and columns SNP (rsID) BETA A1 (effect allele) A2 and P. This will be the input for multivariable LDSR ( Figure 1 ). Taking the files produced in step 4 (ending .gz), use the GenomicSEM R package to carry out multivariable LDSR on all the traits you intend to include in the common-factor GWAS. Multivariable LDSR is a means to calculate genetic correlation between all traits, and produces the empirical genetic covariance matrix used for estimating parameters of genomic SEMs. Note: Be sure to assign the sample and population prevalence values (10 & 11), and to maintain trait order and order of corresponding attributes (sample prevalence, population prevalence, trait names). Example of files for GenomicSEM ‘sumstats’ function and files for input in GenomicSEM or ldsc ‘munge’ Timing: approximately 1 h (parallel, 1 h per chromosome, all chromosomes run at once) to 1 day (sequential) Timing: approximately 1 h (parallel, 1 h per chromosome, all chromosomes run at once) to 1 day (sequential) Here we describe how to apply ‘sumstats’ to summary statistics. This function performs final cleaning and harmonization steps across all included GWAS traits within chromosome. 15. Use the GenomicSEM R package ‘sumstats’ function. The input for this step is sumstats files that you split by chromosome (22 files per trait). Note: This step can in theory be carried out genome-wide (i.e., without first splitting the summary statistic files by chromosome into 22 separate files) but this usually increases computing walltime significantly, and since each chromosome run is independent of all others, it can be efficient to run this step in parallel especially on a personal laptop. Note also that when the sample size is already provided in the summary statistics files, ‘N’ object is set to value NULL rather than a vector of sample sizes. Use the GenomicSEM R package ‘sumstats’ function. The input for this step is sumstats files that you split by chromosome (22 files per trait). Note: This step can in theory be carried out genome-wide (i.e., without first splitting the summary statistic files by chromosome into 22 separate files) but this usually increases computing walltime significantly, and since each chromosome run is independent of all others, it can be efficient to run this step in parallel especially on a personal laptop. Note also that when the sample size is already provided in the summary statistics files, ‘N’ object is set to value NULL rather than a vector of sample sizes. Timing: 1–5 h Timing: 1–5 h Here we describe how to run a common-factor GWAS model without including individual SNP effects. 16. Using the GenomicSEM R package ready-made ‘commonfactor’ function and the multivariable LDSR output (step 5), run the common-factor GWAS model and save the output. 17. Assess the model fit – this can be accessed from the output object with the addition of ‘$modelfit’ and/or ‘$results’. If each trait is significantly associated with the latent common factor and the relationships between traits are broadly as expected, proceed to fitting the model with individual SNP effects (step 18). Take note of any warnings as some contain instructions on options for optimizing the model. Note: To a degree this is subjective and requires consideration, and background understanding of traits and GWAS included. In some cases, you may prefer to first fit a large number of different models and formally compare model fit (i.e., exploratory factor analysis) – examples of exploratory factor analyses and how to judge model fit are discussed in-depth in the original GenomicSEM publication 2 (see in particular Online Supplement pg. 8 ‘Model comparisons: neuroticism example’) . In our case we decided on the model a priori , and found the model fit acceptable based on the fact each trait was significantly associated with the latent common factor. In an exploratory factor analysis context you could fit a range of different models that you hypothesize are good candidates to represent between-trait relationships, represented visually by different path diagrams, and then compare fit indices (chi-squared, Aikake Information Criterion (AIC), Standardized Root Mean-Square Residual (SRMR), Comparative Fit Index (CFI)) of each in order to select a best-fitting model. Using the GenomicSEM R package ready-made ‘commonfactor’ function and the multivariable LDSR output (step 5), run the common-factor GWAS model and save the output. Assess the model fit – this can be accessed from the output object with the addition of ‘$modelfit’ and/or ‘$results’. If each trait is significantly associated with the latent common factor and the relationships between traits are broadly as expected, proceed to fitting the model with individual SNP effects (step 18). Take note of any warnings as some contain instructions on options for optimizing the model. Note: To a degree this is subjective and requires consideration, and background understanding of traits and GWAS included. In some cases, you may prefer to first fit a large number of different models and formally compare model fit (i.e., exploratory factor analysis) – examples of exploratory factor analyses and how to judge model fit are discussed in-depth in the original GenomicSEM publication 2 (see in particular Online Supplement pg. 8 ‘Model comparisons: neuroticism example’) . In our case we decided on the model a priori , and found the model fit acceptable based on the fact each trait was significantly associated with the latent common factor. In an exploratory factor analysis context you could fit a range of different models that you hypothesize are good candidates to represent between-trait relationships, represented visually by different path diagrams, and then compare fit indices (chi-squared, Aikake Information Criterion (AIC), Standardized Root Mean-Square Residual (SRMR), Comparative Fit Index (CFI)) of each in order to select a best-fitting model. Timing: this is dependent on compute resources, model complexity, number of traits, and number of SNPs; estimate from an 8-core laptop, 5 indicators (traits), and 1 million SNPs: 8 h. On our HPC, 28 cores, 6 indicators: 1 h per chromosome Timing: this is dependent on compute resources, model complexity, number of traits, and number of SNPs; estimate from an 8-core laptop, 5 indicators (traits), and 1 million SNPs: 8 h. On our HPC, 28 cores, 6 indicators: 1 h per chromosome Here we describe steps to run the common-factor GWAS (applying the SEM and including individual SNP effects to estimate association between each SNP and the common latent factor). 18. Using GenomicSEM R package ‘commonfactor’ function, multivariable LDSR output, and the output from ‘sumstats’ per chromosome, carry out common-factor GWAS separately per chromosome. Note: The resulting output will be collated afterwards. Again, this can be conducted across all chromosomes at once, but this is time and computation-intensive. The output files will look similar to standard GWAS summary statistics, with 1 row per SNP, a CHR, BP, A1,A2 column, and additional columns for values related to estimates of heterogeneity, whether the model produced warnings for that SNP, and whether the model fails to converge for that SNP. 19. Calculate a sample size to accompany the common-factor GWAS output (Ñ). This estimate is calculated by first restricting the common-factor GWAS output to SNPs with MAF 0.1, then finding, per-SNP: (Equation 3) 1 2 ∗ M A F ∗ ( 1 − M A F ) ∗ S E 2 Using GenomicSEM R package ‘commonfactor’ function, multivariable LDSR output, and the output from ‘sumstats’ per chromosome, carry out common-factor GWAS separately per chromosome. Note: The resulting output will be collated afterwards. Again, this can be conducted across all chromosomes at once, but this is time and computation-intensive. The output files will look similar to standard GWAS summary statistics, with 1 row per SNP, a CHR, BP, A1,A2 column, and additional columns for values related to estimates of heterogeneity, whether the model produced warnings for that SNP, and whether the model fails to converge for that SNP. Calculate a sample size to accompany the common-factor GWAS output (Ñ). This estimate is calculated by first restricting the common-factor GWAS output to SNPs with MAF 0.1, then finding, per-SNP: The mean of these values then gives Ñ. This value can be rounded down to a whole number where necessary (e.g., LD score regression to estimate genetic overlap with other traits requires whole number N values for each input trait GWAS). Note: This will be necessary for any downstream analyses like LDSR, and is useful to know to assess how powerful the common-factor GWAS ended up being. Restricting SNPs included in this calculation to those with MAF > 0.1 and < 0.4 is based on analyses by Mallard et al. 12 (see their Supplement) and produces sample size estimates that are more stable (compared to those produced when SNPs of all MAF ranges are included). Note: This will be necessary for any downstream analyses like LDSR, and is useful to know to assess how powerful the common-factor GWAS ended up being. Restricting SNPs included in this calculation to those with MAF > 0.1 and < 0.4 is based on analyses by Mallard et al. 12 (see their Supplement) and produces sample size estimates that are more stable (compared to those produced when SNPs of all MAF ranges are included).

The authors declare no competing interests.

L.M.H. acknowledges funding from NIMH ( R01MH124839 , R01MH118278 , R01MH125938 , RM1MH132648 , and R01MH136149 ), NIEHS ( R01ES033630 ), and the Department of Defense ( TP220451 ). Core facilities used include the Yale McCleary HPC cluster.

Unusual or unexpected genetic correlations may be observed between traits (related to step 14). For example, we initially observed unexpected negative genetic correlation between endometriosis and all other COPCs in the multivariable LDSR output (step 5). This was due to misinterpretation of the effect vs non-effect allele in the GWAS summary statistics, and corrected after we consulted the FinnGen website. • Always check what the GWAS summary statistic file columns refer to and how they should be interpreted. What is the effect allele and how is it labeled? If the effect allele is mismatched for one GWAS trait compared to the others, ‘swap’ effect alleles to match, and multiply beta (effect size) by −1. Always check what the GWAS summary statistic file columns refer to and how they should be interpreted. What is the effect allele and how is it labeled? If the effect allele is mismatched for one GWAS trait compared to the others, ‘swap’ effect alleles to match, and multiply beta (effect size) by −1. File order must be maintained throughout all steps of the analysis (related to steps 13 onward). For example, the order of munged sumstat files fed to the ‘ldsc’ function for multivariable LDSC must match up with the order of sample prevalence values given to the ‘sample.prev’ object, and this must also match the values in the ‘population.prev’ object, and all this needs to correspond to the correct trait name in ‘trait.names’. Later on, the order of elements in ‘trait.names’, ‘se.logit’, ‘files’,’OLS’, ‘linprob’ must all match, and be in the same order as the order of traits fed to ldsc previously. • Choose an order of traits, note this, and stick to it throughout. Instead of manually running a script per chromosome (if splitting analyses per-chromosome), make use of loops and other methods to methodically and uniformly carry out analyses over chromosomes (e.g., making file lists to use inside other scripts, in the file order decided at the start). chr<-1:22 rm(i) for(i in 1:length(chr)){ c1<-chr[i] f1<-c(sprintf("CWP_GWAS_EU_ANCESTRY_UKB_rsid_maf_chr%s",c1), sprintf("IBS_GCST90016564_buildGRCh37_Neff_chr%s",c1), sprintf("headache2_2017-10-12.allele_labels_chr%s",c1), sprintf("phenotype_1_GWA_metal_result_Geisinger_eMERGE1.TBL.rsids.maf_chr%s",c1), sprintf("TMJ_GCST90044107_buildGRCh37_maf_chr%s", c1), sprintf("finngen_R9_N14_ENDOMETRIOSIS_NOS_betaFix_hg19_preMunge_chr%s", c1)) filename<-sprintf("new_filelist_noMIG_chr%s",c1) write.table(f1, file=filename, row.names=F, col.names=F, quote=F, sep="∖t") rm(c1) rm(f1) rm(filename) } Choose an order of traits, note this, and stick to it throughout. Instead of manually running a script per chromosome (if splitting analyses per-chromosome), make use of loops and other methods to methodically and uniformly carry out analyses over chromosomes (e.g., making file lists to use inside other scripts, in the file order decided at the start). chr<-1:22 rm(i) for(i in 1:length(chr)){ c1<-chr[i] f1<-c(sprintf("CWP_GWAS_EU_ANCESTRY_UKB_rsid_maf_chr%s",c1), sprintf("IBS_GCST90016564_buildGRCh37_Neff_chr%s",c1), sprintf("headache2_2017-10-12.allele_labels_chr%s",c1), sprintf("phenotype_1_GWA_metal_result_Geisinger_eMERGE1.TBL.rsids.maf_chr%s",c1), sprintf("TMJ_GCST90044107_buildGRCh37_maf_chr%s", c1), sprintf("finngen_R9_N14_ENDOMETRIOSIS_NOS_betaFix_hg19_preMunge_chr%s", c1)) filename<-sprintf("new_filelist_noMIG_chr%s",c1) write.table(f1, file=filename, row.names=F, col.names=F, quote=F, sep="∖t") rm(c1) rm(f1) rm(filename) } Common-factor GWAS runs very slowly when queueing R scripts in a Linux environment (related to step 18). This is a known issue and explained further in the GenomicSEM GitHub wiki (1.1). This is caused by instances of R on Linux where a parallel backend is configured to use the maximum number of cores, which when combined with processes which are part of GenomicSEM, generates too many threads and leads to congestion. • This can be fixed using the code recommended in the GenomicSEM wiki, which we specifically used in .sh scheduling scripts. This can be added into whatever shell script is used to queue and run the common-factor GWAS R script, and should be put before any ‘Rscript’ command. export OPENBLAS_NUM_THREADS=1 OMP_NUM_THREADS=1 NUMEXPR_NUM_THREADS=1 VECLIB_MAXIMUM_THREADS=1 This can be fixed using the code recommended in the GenomicSEM wiki, which we specifically used in .sh scheduling scripts. This can be added into whatever shell script is used to queue and run the common-factor GWAS R script, and should be put before any ‘Rscript’ command. export OPENBLAS_NUM_THREADS=1 OMP_NUM_THREADS=1 NUMEXPR_NUM_THREADS=1 VECLIB_MAXIMUM_THREADS=1 There may be various errors introduced during the ‘munge’ and ‘sumstats’ function steps, due to column name ambiguity (related to steps 13-15). If there are multiple columns in the input summary statistic files for munge/sumstats functions with similar names (e.g., N column and an Neff column in the same file), this can lead to confusion, as can non-standard column names (e.g., names like rs_number instead of ‘SNP’ as a column name). • Column names you wish to include in the munge/sumstats functions can be specified as part of the function, alternatively excess unused columns can be removed before summary statistic files are saved (‘prepare summary statistics’). Column names you wish to include in the munge/sumstats functions can be specified as part of the function, alternatively excess unused columns can be removed before summary statistic files are saved (‘prepare summary statistics’). If a well-powered, appropriate GWAS of your trait(s) of interest is not used, common-factor GWAS findings may be misleading (related to steps 3 & 4). • It is good practice to ensure all GWAS in the analysis were carried out on the same genetic ancestry group. • Additionally, ensure that sample size is sufficiently large. Sample sizes needed to achieve sufficient power to find trait-associated SNPs in GWAS have been discussed extensively, 15 , 16 , 17 and will depend on characteristics of the trait(s). These characteristics include whether traits are case/control (binary) or quantitative, their genetic architecture in terms of heritability and degree of polygenicity, and sample prevalence of the binary trait. Study design qualities such as whether case/control ratios balanced and what underlying model has been chosen for the GWAS will also influence power. As a general rough estimate, GWAS sample sizes in well-powered studies tend to range from tens of thousands to millions. • Check that the type of genotyping array is appropriate - for instance, we identified a GWAS of myalgic encephalitis, 18 but this was performed using a specific genotyping array for putative immune variants rather than an assessment of genome-wide genotyped variation, so this was not included. It is good practice to ensure all GWAS in the analysis were carried out on the same genetic ancestry group. Additionally, ensure that sample size is sufficiently large. Sample sizes needed to achieve sufficient power to find trait-associated SNPs in GWAS have been discussed extensively, 15 , 16 , 17 and will depend on characteristics of the trait(s). These characteristics include whether traits are case/control (binary) or quantitative, their genetic architecture in terms of heritability and degree of polygenicity, and sample prevalence of the binary trait. Study design qualities such as whether case/control ratios balanced and what underlying model has been chosen for the GWAS will also influence power. As a general rough estimate, GWAS sample sizes in well-powered studies tend to range from tens of thousands to millions. Check that the type of genotyping array is appropriate - for instance, we identified a GWAS of myalgic encephalitis, 18 but this was performed using a specific genotyping array for putative immune variants rather than an assessment of genome-wide genotyped variation, so this was not included.