{'fraction_non_alphanumeric': 0.03376371650983212, 'fraction_numerical': 0.027417388314002564, 'mean_word_length': 5.083682008368201, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A scalable Bayesian functional GWAS method accounting for multivariate quantitative functional annotations with applications for studying Alzheimer disease

Summary Existing methods for integrating functional annotations in genome-wide association studies (GWASs) to fine-map and prioritize potential causal variants are limited to using non-overlapped categorical annotations or limited by the computation burden of modeling genome-wide variants. To overcome these limitations, we propose a scalable Bayesian functional GWAS method to account for multivariate quantitative functional annotations (BFGWAS_QUANT), accompanied by a scalable computation algorithm enabling joint modeling of genome-wide variants. Simulation studies validated the performance of BFGWAS_QUANT for accurately quantifying annotation enrichment and improving GWAS power. Applying BFGWAS_QUANT to study five Alzheimer disease (AD)-related phenotypes using individual-level GWAS data (n = ∼1,000), we found that histone modification annotations have higher enrichment than expression quantitative trait locus (eQTL) annotations for all considered phenotypes, with the highest enrichment in H3K27me3 (polycomb regression). We also found that cis-eQTLs in microglia had higher enrichment than eQTLs of bulk brain frontal cortex tissue for all considered phenotypes. A similar enrichment pattern was also identified using the International Genomics of Alzheimer’s Project (IGAP) summary-level GWAS data of AD (n = ∼54,000). The strongest known APOE E4 risk allele was identified for all five phenotypes, and the APOE locus was validated using the IGAP data. BFGWAS_QUANT fine-mapped 32 significant variants from 1,073 genome-wide significant variants in the IGAP data. We also demonstrated that the polygenic risk scores (PRSs) using effect size estimates by BFGWAS_QUANT had a similar prediction accuracy as other methods assuming a sparse causal model. Overall, BFGWAS_QUANT is a useful GWAS tool for quantifying annotation enrichment and prioritizing potential causal variants.


Introduction
Although thousands of significant associations have been identified by single-variant genome-wide association studies (GWASs) for complex traits and diseases, the majority of GWAS signals reside in the noncoding genome regions and have unknown biological meaning. [1][2][3] Existing GWAS results based on single-variant tests are still difficult to interpret with respect to the underlying biological mechanisms. 4,5 Promising advancements in sequencing technology have made plenteous multi-omics data available that provide functional annotations of genetic variants available to the scientific community: Combined Annotation-Dependent Depletion (CADD) score; 6 the Roadmap Epigenomics Mapping Consortium, 7 providing DNA methylation, histone modification, and chromatin accessibility information for various human tissues; the Encyclopedia of DNA Elements (ENCODE), 8 providing functional information of human and mouse genomes; and the Genotype-Tissue Expression (GTEx) project, 9 providing expression quantitative trait locus (eQTL) information of 54 human tissues. In particular, molecular QTLs mapped from profiles of molecular phenotypes (e.g., gene expression from GTEx, 9 chromatin marks from Roadmap, 7,10,11 and protein abundances 12 ) and corresponding genomes (genotype data) have been shown to be enriched with GWAS signals and help interpretate the underlying biology for studying complex traits and diseases. [13][14][15] These molecular QTL have been leveraged especially to prioritize GWAS associations of complex traits and diseases. [16][17][18][19][20] An intuitive but widely used ad hoc approach is to finemap and prioritize potential causal GWAS signals that are also molecular QTLs 21 or located in a region with histone modifications. Recently, advanced statistical methods have been proposed to integrate non-overlapped categorical functional annotation (assigning one function label per variant) with GWAS data to fine-map GWAS results. [22][23][24][25][26] PAINTOR has been proposed to integrate multivariate quantitative functional annotations with GWAS summary statistics to fine-map GWAS loci with thousands of variants 27,28 ; FunSPU 29 and STAAR 30 have been proposed to incorporate multiple biological annotations for rare variant association tests. These existing methods have shown the feasibility and promising results of integrating multivariate quantitative functional annotations with GWAS data to fine-map GWAS results and prioritize potential causal variants. However, these methods were not developed for joint modeling millions of genome-wide variants with multivariate quantitative annotations, which would lead to less accurate quantification of annotation enrichment and reduced power of fine-mapping.
Here, we propose a scalable Bayesian functional GWAS method for integrating multivariate quantitative functional annotations with GWAS data by a Bayesian hierarchical variable selection regression model, referred to as BFGWAS_QUANT. BFGWAS_QUANT assumes a hierarchical logistic prior for the causal probabilities of genetic variants in the standard Bayesian variable selection regression (BVSR)-based GWAS method 31 to jointly model millions of genome-wide genetic variants. BFGWAS_QUANT adapts the scalable expectation maximization Markov chain Monte Carlo (EM-MCMC) algorithm developed by the previous Bayesian functional GWAS (BFGWAS) method, which only models non-overlapped categorical annotations (referred to as BFGWAS_CAT in this paper). 22 BFGWAS_QUANT further improves the computation efficiency by pre-calculating the linkage disequilibrium (LD) correlation matrix and single variant test Z score statistics that are used in the MCMC algorithm or using reference LD correlation matrix and summary-level GWAS data. Bayesian causal posterior probability (CPP) and genetic effect size estimates will be generated by BFGWAS_QUANT, along with enrichment quantification of considered multivariate quantitative functional annotations. Bayesian estimates of genetic effect sizes can be used to derive polygenic risk scores (PRSs) that account for functional annotations.
By simulation studies, we showed that our Bayesian estimates of functional enrichment converged and GWAS power was improved over the standard BVSR method without accounting for functional annotations. 31 We then applied BFGWAS_QUANT to real GWAS data for studying Alzheimer disease (AD) related phenotypes, 32,33 accounting for multivariate quantitative annotations with respect to Roadmap histone modifications 7 of brain mid-frontal gyrus, eQTLs of brain frontal cortex tissue 32,34 and eQTLs of microglia cell type. 35,36 We showed that BFGWAS_ QUANT identified interesting enrichment patterns and generated fine-mapped GWAS results using individuallevel and summary-level GWAS data. We also showed that PRSs derived from BFGWAS_QUANT effect size estimates led to similar accurate AD risk prediction as other PRS methods assuming a sparse causal model.
Under Material and methods, we provide an overview of the BFGWAS_QUANT method, Roadmap histone modification and eQTL-based quantitative functional annotations, simulation study design, and application studies of AD. Under Results, we describe the results of simulation and application studies. We then end with a Discussion of the advantages and limitations of BFGWAS_QUANT.

Material and methods
Hierarchical BVSR model BFGWAS_QUANT assumes a hierarchical BVSR 31 model for genome-wide variants, where y n31 is a vector of standardized phenotype with n subjects, X n3p is the standardized genotype matrix with p genome-wide genetic variants, and b p31 is the vector of the genetic effect sizes. Spike-and-slab variable selection prior is assumed per effect size b i . That is, b i has probability p i to be non-zero and follows a normal distribution centered at zero and probability ð1 À p i Þ to be zero with a point-mass density function at 0, where p i denotes the ''casual'' probability of the i th variant. We assume a hierarchical logistic model for the causal probabilities of genetic variants to account for multivariate quantitative functional annotations, where p i denotes the casual probability of the i th variant as in .; A iJ Þ 0 denotes the augmented annotation vector for the i th variant with an intercept term as the first element, and coefficient vector a ¼ ða 0 ; a 1 ; .; a J Þ 0 denotes the intercept term a 0 and enrichment quantification with respect to functional annotation j ¼ 1; .; J. Further, we assume a fixed value in the domain of (0, 1] for t b ; a fixed value for a 0 in the domain of (À13.8, À9), and a standard normal prior for enrichment parameters ða j $ Nð0; 1Þ; j ¼ 0; 1; .JÞ. In particular, t b ¼ 1 would assume that the prior variance of effect sizes is the same as the marginal effect size estimates in a single-variant regression model, and smaller t b values would inflate the magnitude of Bayesian effect size estimates. The lower bound value of a 0 ¼ À 13:8 would assume that the prior causal probability is 10 À6 when a j ¼ 0; j ¼ 1; .; J (see supplemental information for model details).

Adapted EM-MCMC algorithm
To overcome the heavy computational burden and poor mixing rate of posterior samplings by the standard MCMC algorithm, 31 we adapt the scalable EM-MCMC algorithm developed for the original BFGWAS method. 22 Specifically, we first segment genome-wide variants into approximately independent genome blocks with about 5,000-10,000 variants based on LD structure. 22,37 Second, conditioning on given hyperparameters ða j ; j ¼ 1; .; JÞ, we conduct a standard MCMC algorithm within each genome block to obtain Bayesian posterior estimates for genetic effect size (b i ) and causal probability (p i ) per variant (expectation (E) step). Third, conditioning on the Bayesian estimates of ðb i ; p i ; i ¼ 1; .; pÞ, we update the values of hyperparameters ða j ; i ¼ 1; .; JÞ by maximizing their conditional posterior likelihood (maximization [M] step). The computational optimization algorithm Broyden-Fletcher-Goldfarb-Shanno (BFGS) 38 is used to obtain the maximum a posteriori (MAP) estimates for a j ; j ¼ 1; .; J. The EM steps will be iterated ($5 iterations) until the estimates of hyperparameters converge.
The Bayesian estimate for the causal probability (p i ) per variant is referred to as CPP, and the Bayesian estimates for annotation coefficients ða j ; j ¼ 1; .; JÞ are referred to as quantified enrichment of multivariate functional annotations. SNPs with CPP greater than 0.1068 will be considered significantly associated with the phenotype of interest, where the significance threshold has been shown to be equivalent to p < 5 3 10 À8 in a previously published BFGWAS_CAT paper. 22 Because a multivariate regression model is fitted per genome block, the LD among all variants per genome block is accounted for during the Bayesian inference of CPP and effect size b. The GWAS results obtained by BFGWAS_QUANT will be fine-mapped, and variants with enriched annotations will be prioritized.
In particular, implementing the MCMC algorithm per genome block can greatly reduce the search space (from genome-wide to a genome block) and facilitate parallel computing (one core per genome block), leading to an efficient convergence rate and improved mixing rate. Computation efficiency is further improved by implementing the MCMC algorithm using a precalculated LD correlation matrix per genome block and singlevariant test Z score statistics or a reference LD correlation matrix and GWAS Z score statistics, which will save up to 90% computation time compared with using individual-level GWAS data. 37 With 32 computation cores in one node, BFGWAS_QUANT can complete analyzing approximately 10 million SNPs in approximately 4 h for 5 EM iterations.

eQTL-based functional annotations
In this paper, we considered 5 real eQTL-based quantitative functional annotations. Three of these annotations (Allcis-eQTL, 95% CredibleSet, and MaxCPP) were constructed based on standard cis-eQTL data of brain frontal cortex tissue from GTEx data 9,34 : (1) binary annotation Allcis-eQTL was constructed by taking all SNPs that were identified as a significant cis-eQTL (false discovery rate [FDR] < 5%; 1 Mb from the transcription start site [TSS]) for at least one expression quantitative trait (one gene) as 1 and otherwise as 0. (2) For each gene expression trait that has at least one significant (FDR < 5%) cis-eQTL, CAVIAR 16 was used to calculate the CPP (cis-CPP) of each cis-SNP and identify a 95% credible set. SNPs that do not belong to any 95% credible set were taken as 0 or otherwise 1 for the annotation of 95%Cred-ibleSet. (3) Maximum cis-CPPs per SNP across all genes were taken as quantitative values of MaxCPP. We also took the maximum Bayesian genome-wide CPP of being cis-or trans-eQTL across all genes in brain frontal cortex tissue from the Religious Orders Study and Rush Memory and Aging Project (ROS/MAP) 21,33 as the fourth annotation BGW_MaxCPP, where cis-and trans-CPPs were estimated by Bayesian genome-wide transcriptome-wide association study (BGW-TWAS) method. 32,37 Last, we derived a fifth Microglia-eQTL annotation from two datasets of recent microglia cell-type specific eQTL summary statistics, 35,36 where 1 indicates being identified as a cis-eQTL to any gene in either microglia dataset or otherwise 0.

Histone modification-based functional annotation
We constructed 5 histone modification-based functional annotations using the epigenomics data of core histone modifications in the brain mid-frontal gyrus region from the Roadmap Epigenomics database: 7 H3K4me1 (primed enhancers), H3K4me3 (promoters), H3K36me3 (gene bodies), H3K27me3 (polycomb regression), and H3K9me3 (heterochromatin). For each histone modification, peak regions from replicates of the same sample were first merged and then overlapped with peak regions of other samples by Bedtools (v.2.27.0). 39 If a genetic variant resides in the overlapped peak regions of a histone modification, then 1 would be assigned to the function annotation of such a modification for this variant, or 0 would be assigned otherwise.

Simulation study design
We conducted simulation studies to validate the performance of BFGWAS_QUANT. Continuous phenotypes were simulated using the real whole-genome sequence (WGS) data of chromosomes 19 (122,745 SNPs with minor allele frequency [MAF] > 0.01) for 1,893 samples from the ROS/MAP cohort 33,40 and Mount Sinai Brain Bank (MSBB) study. 41 Phenotypes were simulated based on the multivariate linear additive model (Equation 1) with true genetic effect sizes b p31 generated based on the hierarchical logistic model with multivariate quantitative functional annotations (Equation 2). Scenarios with various numbers of true causals and heritability were considered.
Besides the real cis-eQTL-based functional annotations of Allcis-QTL, 95%CredibleSet, and MaxCPP, we also considered a fourth artificial annotation randomly generated from N(0, 1) as a negative control. With chosen annotation enrichment parameters (a 0 ¼ ð À 10:5; À 9:5Þ;a 1 ¼ 4;a 2 ¼ 1:5;a 3 ¼ 0:5;a 4 ¼ 0), we first calculated casual probabilities ðp i Þ for all considered 122,745 SNPs by Equation 2, where a 0 was chosen to ensure the total number of true causal SNPs fall in (5, 10) with a 0 ¼ À 10:5 or (15, 30) with a 0 ¼ À 9:5. Second, a vector of binary indicator ðg i Þ of true causal SNPs was generated from the corresponding Bernoulli distribution with probability (p i ) for i ¼ 1; .; p. Third, genetic effect sizes were taken as 0 for SNPs with g i ¼ 0 or generated from a normal distribution for SNPs with g i ¼ 1. Finally, phenotypes were generated from Equation 1 with simulated genetic effect sizes and random errors e $ Nð0; ð1 À h 2 ÞIÞ to ensure that a target total heritability h 2 ¼ ð0:25; 0:5Þ was equally explained by all true causal SNPs. Four scenarios were considered, including one with relatively sparse true causals in the range of (5, 10) and one with the number of true causals in the range of (15, 30) with respect to two different heritability values (0.25, 0.5).
We considered a null enrichment scenario where none of the annotations were enriched. In this scenario, we randomly selected 10 true causal SNPs and assigned them genetic effect sizes generated from a normal distribution. We simulated the phenotype as described above with a targeted h 2 ¼ 0.5.
We repeated 100 simulations per scenario to evaluate our Bayesian estimates for annotation enrichment, total heritability, and true causal SNPs with respect to sensitivity (power) and positive predictive values (PPVs). 42 Sensitivity (power) is defined as the proportion of true positive findings among all true casual variants, and PPV is defined as the proportion of true positive findings among all identified significant associations. We took simulated true casual SNPs and those having R 2 > 0:3 with true casual SNPs with b p i >0.1 as true positive findings, following the significance rule used by the BFGWAS_CAT method. 22 The sensitivity and PPV are given by We compared this with the standard BVSR method, 31 which does not account for functional annotations. We estimated the total heritability by the squared correlation between the simulated phenotypes and the PRSs based on Bayesian estimates of genetic Applications to ROS/MAP individual-level GWAS data ROS/MAP are two prospective longitudinal community-based cohort studies that recruit older adults without known dementia at baseline and follow up with participants annually until the time of death. 33,40,43 Participants agree to annual clinical evaluation and brain autopsy at the time of death, signing an informed consent form and Anatomic Gift Act. All participants in this study also sign a repository consent to allow their data to be re-purposed. WGS data were profiled for 1,200 samples by using the KAPA Hyper Library Preparation Kit and Illumina HiSeq X sequencer.
We applied BFGWAS_QUANT to account for five eQTL-based and five histone modification-based quantitative annotations as described above to study five AD-related phenotypes, 9,[32][33][34] including the binary clinical diagnosis of late-onset Alzheimer dementia (n ¼ 1,087), three quantified postmortem pathology indices of AD (i.e., PHFtau tangle density with n ¼ 1,105, b-amyloid load with n ¼ 1,113, and a global measurement of AD pathol-ogy burden with n ¼ 1,123), and a quantitative measurement of cognition decline rate with n ¼ 1,049. 33,40 The cognition decline rate was constructed as the random slope per sample from a linear mixed model of annual longitudinal measurements of cognition function. Details about ROS/MAP and phenotypic traits measured can be found in a previously published paper. 44 We also adjusted for the covariates of age, sex, smoking status, study index (ROS or MAP), and first 3 genotype principal components by regressing these covariates out from the phenotypes and taking the corresponding regression residuals as the outcome in the BFGWAS_ QUANT method.

Application to IGAP summary-level GWAS data
We applied BFGWAS_QUANT to the stage 1 summary-level GWAS data of the International Genomics of Alzheimer's Project (IGAP), 45 along with the above 10 eQTL-and histone modification-based functional annotations and reference LD generated from ROS/MAP. The stage 1 IGAP summary-level GWAS data were generated by meta-analyses consisting of 17,008 individuals with AD and 37,154 control individuals of European ancestry (n ¼ $54,000).

AD risk prediction by PRS
To show the usefulness of PRSs for risk prediction, we evaluated two sets of PRSs that were, respectively, derived from BFGWAS_QUANT and BVSR summary statistics using ROS/MAP GWAS data. We used the independent test samples from Mayo Clinic Alzheimer's Disease Genetics Studies (MCADGS). 46,47 MCADGS contain 2,099 European-descent samples (844 individuals with AD and 1,255 control individuals) with microarray genotype data profiled that were further imputed to the 1000 Genome Project Phase. 48 We compared PRSs using Bayesian effect size estimates with three commonly used PRS methods: the standard method using informed LD pruning and p value thresholding (P þ T), 49,50 LDpred2, 51 and PRS-CS. 52 p value thresholds of ð10 À 2 ; 10 À 3 ; 10 À 5 ; 5 310 À 8 ; 10 À 8 Þ and LD thresholds of ð0:1; 0:3; 0:5; 0:7; 0:9Þ were considered by the P þ T method. Reference LD derived from 1000 Genome data were used by the LDpred and PRS-CS methods. AD risk prediction accuracy was evaluated using the area under the receiver operating characteristic (ROC) 53 curve (AUC) for MCADGC test samples.

Ethics statement
The ROS/MAP and MCADGS data analyzed in this study were generated with approval of the institutional review board (IRB) of Rush University Medical Center, Chicago, IL, and Mayo Clinic, respectively. All samples analyzed in this study were de-identified, and all analyses were approved by the IRB of Emory University School of Medicine.

Simulation results
For all considered simulation scenarios, our Bayesian estimates of functional annotation enrichment achieved convergence with 4 EM iterations, as shown in boxplots of 100 simulation replicates (Figures 1A, 1D, and S1). Although BFGWAS_QUANT overestimated a 1 and a 2 , the Bayesian estimates still reflect the correct enrichment pattern among all considered annotations. By taking the 2.5th and 97.5th quantiles of these 100 Bayesian estimates to estimate the corresponding 2.5th and 97.5th quantiles of the estimator distributions, the true enrichment values ða 1 ¼ 4; a 2 ¼ 1:5; a 3 ¼ 0:5; a 4 ¼ 0Þ indeed fell within this range. For example, in the scenario with a 0 ¼ À 10:5 and h 2 ¼ 0:25 ( Figure 1A), the estimated 2.5th and 97.5th quantiles are ð1:62; 5:50Þ for a 1 , (0.98, 4.57) for a 2 , (0.02, 1.20) for a 3 , and (0.00, 0.38) for a 4 . We observed precise estimates of 0 enrichment for the artificial annotation (Figures 1 and S1) and all annotations in the scenario with null enrichment (Figure S3A), which demonstrated the ability of BFGWAS_QUANT to identify null enrichment.
By taking PRS as estimated phenotypes and taking the squared correlation between RPS and simulated phenotypes as the estimate of phenotype heritability, we obtained similar heritability estimates by BFGWAS_QUANT and BVSR, which are close to the true heritability ( Figures 1B, 1E, and S3B). For scenarios with true enrichment, BFGWAS_QUANT obtained substantially higher sensitivity (power) and similar PPVs compared with BVSR for all scenarios (Figures 1C, 1F, S2, S4, and S5). For the scenario with null enrichment, BFGWAS_QUANT and BVSR performed comparably (Figures S3B-S3D).
These simulation studies validated the usefulness of BFGWAS_QUANT for quantifying multivariate functional annotations, estimating phenotype heritability, and identifying true causal SNPs. By accounting for multivariate quantitative annotations, BFGWAS_QUANT showed improved performance than the standard BVSR method, especially with higher power for identifying true causal SNPs and accurate enrichment estimation.
Application GWAS results for studying AD Applying BFGWAS_QUANT to the individual-level ROS/ MAP and summary-level IGAP GWAS data, we obtained consistent patterns for Bayesian enrichment estimates of 5 eQTL-based and 5 histone modification-based functional annotations (Figures 2, S6, and S7). In particular, the histone modification-based functional annotations had higher enrichment than eQTL-based annotations when studying the individual-level ROS/MAP GWAS data, with the highest enrichment for H3K27me3 and second highest for H3K4me1. BFGWAS_QUANT estimated the second highest enrichment for the Microglia_eQTL annotation when studying the summary-level IGAP GWAS data with a larger sample size (n ¼ $54,000). Even with a small sample size in the individual-level ROS/MAP GWAS data, the Microglia_eQTL annotation was still identified with higher enrichment than other annotations based on eQTL of the bulk brain frontal cortex tissue. There is mounting evidence showing that microglia (composing <10% cells in the bulk brain frontal cortex tissue) play important roles in development and progression of AD pathology, 54 and cell-type-specific differential expression of GWAS risk genes of AD is only present in microglia. 55 Using the ROS/MAP individual-level GWAS data, four significant SNPs with Bayesian CPP greater than 0.1068 were identified for AD (rs429358, rs10414043, rs769449, and rs7256200) by BFGWAS_QUANT (Table 1; Figures 3A, S8, and S9). In particular, SNP rs429358 (CPP ¼ 0.144, p ¼ 7.72 3 10 À13 , missense variant) is the famous known APOE E4 risk allele of AD 56 and has a significant Bayesian CPP greater than 0.1068 for all 5 AD related phenotypes. SNPs rs10414043 (CPP ¼ 0.111, p ¼ 2.71 3 10 À12 ) and rs7256200 (CPP ¼ 0.315, p ¼ 2.71 3 10 À12 , regulatory variant) are upstream of the known risk gene APOC1 of AD and blood protein traits. 37,57,58 Besides the missense APOE E4 risk alleles rs429358 and rs7412, one additional signifi-cant SNP, rs1065853, is identified upstream of APOC1 for global AD pathology, which is a known GWAS signal for blood protein traits such as low-density lipoprotein. 58,59 Of 11 significant SNPs identified for at least one ADrelated phenotype (Table 1), 2 are intergenic, and the other 9 SNPs are intron, regulatory, missense, and stop-gained variants.
Using the IGAP summary-level GWAS data, BFGWAS_ QUANT fine-mapped 32 significant SNPs with Bayesian CPP greater than 0.1068 associated with AD (Table 2; Figure 3B). Multiple SNPs located in genes PVRL2, APOE/ TOMM40, and APOC1 in chromosome 19 were found to be associated with AD with CPP ¼ 1. Interestingly, 10 of 32 significant SNPs are cis-eQTL of microglia, and all significant SNPs except one intergenic SNP (rs7110631) are intron, regulatory, downstream, upstream, 3 0 UTR, and missense variants ( Table 2). All significant SNPs except rs78959900 are located in the peak regions of histone modification H3K27me3 (polycomb regression), which has the highest enrichment. Several SNPs that did not pass the genome-wide significance threshold (p < 5 3 10 À8 ) were identified by integrating 10 functional annotations in BFGWAS_QUANT. These SNPs are located in genes that have been found previously to be genetically linked to AD, such as HLA-DRB1, NUP160, SLC24A4, and CD33. 45,56,[60][61][62][63]  The summation of genome-wide Bayesian CPP of SNPs with CPP greater than 0.01 can be used to estimate the number of total causal SNPs for the phenotype of interest. The threshold of CPP >0.01 is used to exclude adding CPP from random MCMC selections. Although the power is limited for analyzing the ROS/MAP individual-level GWAS data with a small sample size, BFGWAS_QUANT estimated a total of 54 potential causal SNPs for AD using the IGAP summarylevel GWAS data (Table 3).

AD risk prediction by PRS in MCADGS
To show the usefulness of BFGWAS_QUANT for studying complex traits and diseases, we derived a PRS using the Bayesian effect size estimates by BFGWAS_QUANT for an independent GWAS cohort, MCADGS (n ¼ 2,099), and compared the risk prediction accuracy with the PRS using effect size estimates by BVSR, P þ T, LDpred2, and PRS-CS. When the ROS/MAP individual-level GWAS data were used as training data, a comparable AUC was obtained by using Bayesian effect size estimates by BFGWAS_QUANT (0.69) and BVSR (0.68), which was similar to the one obtained by LDpred2-auto (0.68) but significantly higher than the ones obtained by P þ T (0.55), LDpred2-inf (0.53), and PRS-CS (0.54) ( Figure 4A). This showed the advantage of deriving PRSs using Bayesian effect size estimates by BFGWAS_QUANT when individual-level GWAS data are available, especially when the sample size is small. When the IGAP summary-level GWAS data were used as training data, the AUC by BFGWAS_QUANT (0.75) was the same as by LDpred2-auto (0.75) but much lower than P þ T (0.88 with a p value threshold of 10 À 3 ), LDpred2inf (0.94), and PRS-CS (0.93) ( Figure 4B). These results show that an infinitesimal model 64 as assumed by PRS-CS and LDpred2-inf is more suitable for PRS development than the sparse model assumed by BFGWAS_QUANT and LDpred2-auto.

Discussion
We developed a scalable BFGWAS method to account for multivariate quantitative functional annotations (BFGWAS_QUANT) for studying complex traits and diseases based on a hierarchical BVSR model and accompanied by a scalable EM-MCMC computation algorithm. BFGWAS_QUANT has the advantages of quantifying enrichment of functional annotations as well as modeling LD to generate fine-mapped GWAS results that are also prioritized based on their functional annotations. BFGWAS_QUANT can be applied to individual-level and summary-level GWAS data. In particular, the Bayesian effect size estimates can be used to derive a PRS that accounts for functional annotations.
Our simulation studies validated the performance of BFGWAS_QUANT with respect to annotation enrichment quantification, GWAS association identification, and heritability estimation. Compared with BVSR, the BFGWAS_QUANT method had higher sensitivity (i.e., power) and comparable PPV and accuracy of phenotype heritability estimation.
In real studies of AD-related phenotypes using ROS/MAP individual-level and IGAP summary-level GWAS data, we showed that interesting enrichment patterns were identified, fine-mapped GWAS signals were identified, and predictive PRSs were derived. In particular, we found that the histone modification H3K27me3 (polycomb regression) and microglia cis-eQTL annotations were most en-riched for association signals of AD. We also showed that SNPs with single variant test p values < 5 3 10 À 8 could be identified for being prioritized because of their functional annotations.
Despite these advantages, BFGWAS_QUANT does have its limitations. First, the BFGWAS_QUANT model was developed for quantitative traits. However, following previous studies, 22,65 GWAS analysis can still be done for dichotomous traits by quantifying cases as 1 and controls as 0, which will have a similar performance as a probit model when samples are independent and population structure can be addressed by top genotype principal components. Extending the BFGWAS_QUANT method for studying dichotomous traits is also part of our ongoing research. Second, BFGWAS_QUANT assumes the summary-level GWAS data, and reference LDs were derived from populations of the same ancestry. BFGWAS_QUANT per ancestry needs to be applied first, and then the results have to be meta-analyzed for studying GWAS cohorts with  PRSs derived using IGAP summary-level GWAS data . The PRS derived using Bayesian effect size estimates by BFGWAS_QUANT has comparable prediction accuracy as the PRSs derived by BVSR and LDpred2 auto, for all assuming a sparse causal model. PRSs derived by PRS-CS and LDpred2-inf using IGAP summary-level GWAS data as training data have the highest prediction accuracy for assuming an infinitesimal causal model. multiple ancestries. Third, BFGWAS_QUANT assumes a sparse causal genetic architecture that is suitable for generating fine-mapped GWAS results but might lack of power for deriving PRS for complex traits and diseases. Our work demonstrated the usefulness of integrating multivariate quantitative functional annotations in GWASs for quantifying the enrichment of multiple functional annotations and generating fine-mapped GWAS results with higher power. Specifically, accurate quantification of annotation enrichment would help prioritize GWAS signals (fine-mapping) and then help illustrate the underlying genomic etiology of complex traits and diseases. Because publicly available molecular QTL datasets, epigenomic features, and GWAS summary data continue to grow, BFGWAS_QUANT provides a convenient tool for integrative multi-omics analyses of these datasets.