Computation and resource efficient genome-wide association analysis for large-scale imaging studies.

Skip to main content An official website of the United States government Here's how you know Here's how you know Official websites use .gov A .gov website belongs to an official government organization in the United States. Secure .gov websites use HTTPS A lock ( Lock Locked padlock icon ) or https:// means you've safely connected to the .gov website. Share sensitive information only on official, secure websites. Search Log in Dashboard Publications Account settings Log out Search… Search NCBI Primary site navigation Search Logged in as: Dashboard Publications Account settings Log in Search PMC Full-Text Archive Search in PMC Journal List User Guide PERMALINK Copy As a library, NLM provides access to scientific literature. Inclusion in an NLM database does not imply endorsement of, or agreement with, the contents by NLM or the National Institutes of Health. Learn more: PMC Disclaimer | PMC Copyright Notice Nat Commun . 2026 Feb 28;17:3313. doi: 10.1038/s41467-026-69816-z Search in PMC Search in PubMed View in NLM Catalog Add to search Computation and resource efficient genome-wide association analysis for large-scale imaging studies Zhiwen Jiang Zhiwen Jiang 1 Department of Biostatistics, University of North Carolina at Chapel Hill, Chapel Hill, NC USA Find articles by Zhiwen Jiang 1 , Jason Stein Jason Stein 2 Department of Genetics, University of North Carolina at Chapel Hill, Chapel Hill, NC USA Find articles by Jason Stein 2 , Tengfei Li Tengfei Li 3 Department of Radiology, University of North Carolina at Chapel Hill, Chapel Hill, NC USA 4 Biomedical Research Imaging Center, University of North Carolina at Chapel Hill, Chapel Hill, NC USA Find articles by Tengfei Li 3, 4 , Ethan Fang Ethan Fang 5 Department of Biostatistics and Bioinformatics, Duke University, Durham, NC USA Find articles by Ethan Fang 5 , Yun Li Yun Li 1 Department of Biostatistics, University of North Carolina at Chapel Hill, Chapel Hill, NC USA 2 Department of Genetics, University of North Carolina at Chapel Hill, Chapel Hill, NC USA 6 Department of Computer Science, University of North Carolina at Chapel Hill, Chapel Hill, NC USA Find articles by Yun Li 1, 2, 6 , Patrick Sullivan Patrick Sullivan 2 Department of Genetics, University of North Carolina at Chapel Hill, Chapel Hill, NC USA Find articles by Patrick Sullivan 2 , Hongtu Zhu Hongtu Zhu 1 Department of Biostatistics, University of North Carolina at Chapel Hill, Chapel Hill, NC USA 2 Department of Genetics, University of North Carolina at Chapel Hill, Chapel Hill, NC USA 3 Department of Radiology, University of North Carolina at Chapel Hill, Chapel Hill, NC USA 4 Biomedical Research Imaging Center, University of North Carolina at Chapel Hill, Chapel Hill, NC USA 6 Department of Computer Science, University of North Carolina at Chapel Hill, Chapel Hill, NC USA 7 Department of Statistics and Operations Research, University of North Carolina at Chapel Hill, Chapel Hill, NC USA Find articles by Hongtu Zhu 1, 2, 3, 4, 6, 7, ✉ Author information Article notes Copyright and License information 1 Department of Biostatistics, University of North Carolina at Chapel Hill, Chapel Hill, NC USA 2 Department of Genetics, University of North Carolina at Chapel Hill, Chapel Hill, NC USA 3 Department of Radiology, University of North Carolina at Chapel Hill, Chapel Hill, NC USA 4 Biomedical Research Imaging Center, University of North Carolina at Chapel Hill, Chapel Hill, NC USA 5 Department of Biostatistics and Bioinformatics, Duke University, Durham, NC USA 6 Department of Computer Science, University of North Carolina at Chapel Hill, Chapel Hill, NC USA 7 Department of Statistics and Operations Research, University of North Carolina at Chapel Hill, Chapel Hill, NC USA ✉ Corresponding author. Received 2025 Sep 5; Accepted 2026 Feb 9; Collection date 2026. © The Author(s) 2026 Open Access This article is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License, which permits any non-commercial use, sharing, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if you modified the licensed material. You do not have permission under this licence to share adapted material derived from this article or parts of it. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by-nc-nd/4.0/ . PMC Copyright notice PMCID: PMC13066023  PMID: 41764180 Previous version available: This article is based on a previously available preprint posted on medRxiv on November 13, 2025: " Computation and resource efficient genome-wide association analysis for large-scale imaging studies ". Abstract Imaging genetics links genetic variations to brain structures and functions, but the computational challenges posed by high-dimensional imaging and genetic data are significant. In voxel-level genome-wide association studies, we introduce a Representation learning-based Voxel-level Genetic Analysis (RVGA) framework that reduces computational time and storage burden by over 200 times. RVGA enhances statistical power by denoising images and shares minimal datasets of summary statistics for associations across the whole genome of the entire image for secondary analyses. Additionally, it introduces a unified estimator for voxel heritability, genetic correlations between voxels, and cross-trait genetic correlations between voxels and non-imaging phenotypes. Applying RVGA to hippocampus shape and white matter microstructure in the UK Biobank ( n = 53,454) reveals 39 and 275 novel loci, respectively. We identify heterogeneity in heritability within images and subregions that share genetic bases with 14 brain-related phenotypes, such as the genetic correlation between the hippocampus and educational attainment, and between the anterior corona radiata and schizophrenia. RVGA replicates known genetic associations and uncovers new discoveries. Subject terms: Genome-wide association studies, Learning algorithms, Image processing Imaging genetics links genetic variations to brain structures and functions, however the computational challenges posed by high-dimensional imaging and genetic data are significant. Here, the authors develop a scalable genome-wide association analysis framework for imaging data. Introduction Imaging genetics elucidates how genetic variations influence the human brain and links these variations to brain-related disorders such as Alzheimer’s disease 1 – 7 . This field enhances our understanding of the pathophysiological pathways underlying many brain disorders, potentially leading to more precise and personalized treatments. Modern neuroimaging technologies, such as functional magnetic resonance imaging (fMRI), enable us to capture the human brain at high resolution, with 10 4 –10 6 measures per subject. The measure unit can be either voxel (3D volume) or vertex (surface). For simplicity, we use “voxel" to represent both units where no ambiguity arises. However, most existing imaging genome-wide association studies (GWAS) do not fully exploit the rich voxel-level imaging signals and are limited to a few hundred image-derived phenotypes (IDPs), which are generated by aggregating local signals, potentially discarding significant local variations 4 , 8 – 14 . Recent studies have begun to address this issue by employing representation learning to capture essential information from the entire image 15 – 17 and applying the resulting low-dimensional representations (LDRs) in GWAS. However, challenges remain in interpreting the significant loci associated with LDR, reconstructing voxel-level results, and sharing GWAS summary statistics for a large number of IDPs. Voxel-level GWAS (VGWAS) can potentially solve all the limitations, but it faces significant computationally challenges 18 , 19 . First, high-resolution VGWAS generates vast amounts of data, requiring immense computational power and storage capacity, which leads to long processing times and substantial resource demands. Second, sharing VGWAS summary statistics at high resolution is particularly difficult, as their size can be 100 times larger than the NHGRI-EBI GWAS catalog, which contains over 70,000 GWAS datasets. This limits the practical utility and broad dissemination of informative datasets for downstream analyses. Third, correcting for multiple comparisons in hypothesis testing across highly correlated voxels is challenging 20 – 22 . Fourth, noise in imaging data can obscure meaningful genetic associations, reducing statistical power and reliability. This challenge is compounded by the need to accurately identify the genetic architectures underlying complex brain functions and structures, which requires sophisticated statistical methods and large-scale data integration. In contrast, current VGWAS methods, such as considering multiple genetic markers simultaneously 23 , screening candidate markers 24 , 25 , applying rank-reduction 26 – 28 and Bayesian techniques 29 – 31 , or combining multiple test statistics 32 , often overlook the whole-genome summary statistics necessary for secondary analyses, such as heritability and genetic correlations. We introduce the Representation learning-based Voxel-level Genetic Analysis (RVGA) framework to address the aforementioned challenges (Fig. 1 ). First, RVGA enhances both computational and statistical efficiency in VGWAS by decomposing raw images into smooth images and purely random noise, and constructing LDRs to capture essential imaging signals. Second, it proposes a storage-efficient system for sharing summary statistics, enabling all voxel-level secondary analyses. Third, RVGA introduces a unified estimator for voxel heritability, genetic correlations between voxels, and cross-trait genetic correlations between voxels and other phenotypes, using the shared data and a population-matched linkage disequilibrium (LD) matrix. Fig. 1. An overview of RVGA. Open in a new tab A LDRs and functional bases are derived from imaging data using FPCA. The variance-covariance matrix of LDR is estimated by sample covariance. The number of LDRs ( r ) is 1-3 orders of magnitude less than the number of voxels ( N ) because of strong spatial correlation, resulting in substantial dimension reduction. B RVGA conducts VGWAS in two steps. First, LDRs are treated as traits to perform GWAS using external software. Next, VGWAS results are recovered by using the triplets: summary statistics of LDRs, the bases, and the variance-covariance matrix of LDRs. Suppose d SNPs and n subjects in GWAS. The complexity is O ( d n r + d N r ) in time and O ( d r ) in space, which is drastically less than O ( d N n ) and O ( d N ) for scanning each voxel separately. Five advantages identified genetic loci for antages of RVGA over traditional approaches are highlighted. C Using the shared minimal dataset consisting of the triplets, as well as an LD matrix (and its inverse) estimated from a reference panel, RVGA generates an atlas of imaging heritability, estimates genetic correlations between ROIs, and investigates shared genetic bases between images and complex disorders through cross-trait genetic correlation analysis. We used icons from BioRender https://BioRender.com/pd1vhqn . We analyzed hippocampus shape and white matter (WM) microstructure from 53,454 unrelated subjects of European ancestry in the UK Biobank (UKB), splitting into discovery ( n = 33,324) and replication ( n = 20,130). By performing GWAS on 7.8 million common variants across 30,000 hippocampus vertices and 32,217 WM voxels, we identified 39 and 275 novel loci, respectively. We reduced the total file size of summary statistics by 229 times (from 12,611 GB to 55 GB) and shared the data with community. Additionally, we generated atlases of heritability, genetic correlations within regions of interest (ROIs), and cross-trait genetic correlations between ROIs and 14 brain-related phenotypes. Novel genetic correlations were identified between WM tracts and schizophrenia, depression, bipolar disorder, and major depressive disorder. In summary, RVGA surpasses traditional voxel-level imaging genetic frameworks in computation and resource efficiency, facilitates the sharing of VGWAS summary statistics, and pinpoints genetic influences on the human brain and other organs. Results Overview of methods We introduce RVGA, a framework to perform VGWAS and downstream analyses on large-scale imaging genetic datasets. RVGA uses functional principal component analysis (FPCA) to separate smooth imaging signals from white noise and to derive LDRs from images 33 , 34 . These LDRs are used in GWAS to produce whole-genome summary statistics, from which RVGA reconstructs voxel-level summary statistics, creating a comprehensive atlas of genetic associations without the need to store voxel-level data. RVGA offers substantial time and storage savings by 1–3 orders of magnitude. This efficiency not only facilitates data sharing but also enables detailed examination of association patterns across ROIs, enhancing the investigation of genetic architecture at very high resolution. Suppose we have biomedical images for n genetically unrelated subjects (no common ancestor for three generations). Let Y i ( v ) = X i ( v ) + ϵ i ( v ) represent the imaging data for subject i , with N common voxels indexed by v , where X i ( v ) are the underlying images and ϵ i ( v ) are independent white noise terms. The covariance function between two voxels, u and v , C ( u , v ) = Cov( X i ( u ), X i ( v )), has a spectral decomposition C ( u , v ) = ∑ j = 1 ∞ λ j Φ j ( u ) Φ j ( v ) , where λ 1 ⩾ λ 2 ⩾ ⋯ ⩾ 0 are the eigenvalues, and ϕ 1 , ϕ 2 , … are the corresponding functional bases. To estimate C , we first smooth each raw image Y i (⋅) using a local linear estimator 35 with a Gaussian kernel. Each voxel is represented by a weighted average of itself and its neighboring voxels (Methods, Supplementary Note ). We then perform spectral decomposition on the sample covariance of the smoothed images to obtain eigen pairs { λ ^ j , Φ ^ j } j = 1 ∞ . By extracting the top r bases, ϕ ^ r = ( Φ ^ 1 , … , Φ ^ r ) , where r = arg min k ∑ j = 1 k λ ^ j ∑ j = 1 ∞ λ ^ j ≥ q and q represents a proportion of variance (⩾80%), we project the raw image Y i (⋅) onto the functional subspace spanned by ϕ ^ r , yielding LDRs as ξ ^ i j = ⟨ Y i , Φ ^ j ⟩ . The underlying image can be approximated as X i ( v ) ≈ ∑ j = 1 r ξ ^ i j Φ ^ j ( v ) . Thus, the dimensionality of all imaging data Y is reduced from n × N to n × r , with the LDRs Ξ ^ r = ( ξ ^ 1 , … , ξ ^ r ) , as shown in Fig. 1 A. Due to the strong spatial correlation within images, the eigenvalues decay rapidly, making r 1–3 orders of magnitude smaller than N , while still preserving significant information from the original images (Supplementary Fig. 1 ). In GWAS, the marginal genetic effect of a continuous trait is estimated by projecting the trait vector onto the space spanned by the single nucleotide polymorphism (SNP) vector, assuming that all covariate effects have been removed. Similarly, in RVGA, the LDRs are treated as traits for GWAS analysis, projected onto each SNP vector, and the bases Φ ^ r are used to recover the summary statistics for all voxels (Methods, Fig. 1 B). With d SNPs, the computational burden of RVGA is only O ( d n r ) for LDR-level GWAS and O ( d N r ) for recovering VGWAS summary statistics. In contrast, directly performing VGWAS is computationally inefficient, with a complexity of O ( d N n ). In addition to improving computational efficiency, decomposing raw images into underlying smooth images and white noise via the FPCA framework enhances statistical power in GWAS. Importantly, sharing the triplets—the summary statistics of LDRs, the functional bases, and the variance-covariance matrix of LDRs—is sufficient for all secondary analyses (Fig. 1 C). This reduces the storage requirement from O ( d N ) for all VGWAS summary statistics to O ( d r + N r + r 2 ). Next, we introduce a unified estimator to compute voxel SNP heritability, genetic correlations between voxels, and cross-trait genetic correlations with non-imaging phenotypes. Unlike regular definition of heritability, which is the ratio of genetic variance to phenotypic variance, Var( Y (⋅)), voxel heritability is defined as the ratio of genetic variance to the underlying image variance, Var( X (⋅)) (Methods). The estimator is accurate regardless of noise levels. Moreover, with variance analytically derived, it requires only the shared triplets and a population-matched LD matrix (and its inverse). It is far more computationally and statistically efficient than current methods and does not impose assumptions about the distribution or genetic architecture of the effects, making it robust to model mis-specification. Simulation studies Our simulation studies incorporated factors reflecting real imaging genetic studies: potential confounding effects in images, a polynomial decay rate of image eigenvalues (1.8 for slow decay, 2.5 for fast decay), with or without white noise (noise percentage = 50%, 20%, 0%), white noise following Gaussian or Rayleigh distribution, voxel heritability (3%, 10%, 30%), polygenicity (1%, 20% of causal variants), and data truncation levels (75%, 80%, up to 95% of variance preserved by the top LDRs). Refer to the Online Methods section for details about data generation. The relationship between number of LDRs and amount of preserved variance We observed that the amount of variance preserved by the LDRs depends on the decay rate of eigenvalues and noise percentage (Supplementary Fig. 2 ). A slower decay rate indicates reduced spatial correlation, necessitating more LDRs. Noise-free images containing 100% of real imaging signals require more LDRs to capture local variations. Additionally, the number of LDRs remained stable across varying heritability and polygenicity. The accuracy of RVGA GWAS summary statistics We evaluated the accuracy of RVGA GWAS summary statistics compared to VGWAS, using the root-mean-squared-error (RMSE) of z-scores across all voxel-variant pairs. For noise-free images, RMSE of z-scores decreased almost linearly with an increasing number of LDRs (Supplementary Fig. 3 ). For noisy images with Gaussian noise, the most significant RMSE improvement occurred when transitioning from preserving 95% of variance to using raw data (RMSE = 0), and the improvement increased with higher noise level. However, the pattern of RMSE from 95% to 75% of variance was similar across various levels of noise, while noise increased RMSE at all truncation levels. For noisy images with Rayleigh noise, the pattern was similar to the Gaussian case (Supplementary Fig. 4 ). Higher-heritability cases experienced larger bias, but the bias was generally unrelated to polygenicity. Taken together, noise can be removed by using the top LDRs, capturing 95% of variance. RMSE was robust across various proportions of variance, which means in practice, one can flexibly reduce LDRs if decay rate is slow and increase LDRs if computation burden is feasible. The type I error rate and statistical power We evaluated type I error by the proportion of null tests with p -value less than 10 −2 , 10 −3 , and 10 −4 36 . It was well controlled for RVGA regardless of noise distribution, noise percentage, eigenvalue decay rate, and the proportion of variance preserved (Supplementary Figs. 5 – 7 ). We quantified statistical power by the proportion of causal variants identified at the significance level of 10 −4 . There was substantial power gain compared to VGWAS for noisy images (Supplementary Figs. 8 and 9 ). The power for RVGA was unaffected by noise levels and eigenvalue decay rate, while VGWAS was sensitive to noise. That said, RVGA can effectively remove noise accounting for up to 50% of the total variance, across different distributions (normal and Rayleigh), and still achieve comparable power to the noise-free case. indicating RVGA retained image signal. However, we emphasize that selecting only the top PCs cannot preserve all the signal, but the signal loss is minor compared to the noise reduction. As expected, higher heritability and smaller polygenicity (i.e., larger effect size) yielded higher power. We observed slight power inflation in noise-free images with a slow eigenvalue decay rate, which is attributed to the relative smoothness between genetic and non-genetic effects. The heritability and genetic correlation estimator To investigate the robustness of the estimator to various genetic architectures, we varied the coupling strength between genetic effect sizes and minor allele frequency (MAF)/LD using the LDAK model 37 . We also considered different LD matrix types, using either imputed genotype data or genotype array data, with multiple regularization levels. For example, {90%, 85%} denotes preserving 90% of the variance in the LD matrix and 85% in its inverse, achieved by performing eigen-decomposition on each LD block. The heritability estimator was robust to data truncation and polygenicity but sensitive to LD regularization under certain genetic architectures. Specifically, when effect sizes were strongly correlated with both MAF and LD (Supplementary Figs. 10 B and 11 B), restrictive LD regularization led to large mean absolute error (MAE), especially in high-heritability scenarios. This effect was more pronounced in imputed genotype data due to more complex LD between variants. The estimator remained stable in other cases. Using {98%, 95%} for imputed genotype data (Supplementary Figs. 10 and 12 ) and {85%, 80%} for genotype array data (Supplementary Figs. 11 and 13 ) produced unbiased estimates with low MAE, and is thus recommended for real data analysis. We further assessed the impact of white noise on heritability estimation. Applying the estimator to each voxel without noise removal, we observed that noise-laden images resulted in significantly downward-biased estimates, while noise-free images maintained unbiased estimations (Supplementary Figs. 14 and 15 ). However, the estimator was unbiased with noise removal as previously shown. This highlights that RVGA can preserve additive genetic effects while denoising images. Genetic correlation estimates between voxels showed distinct patterns. Data truncation level and true heritability were the main factors affecting MAE, while polygenicity and genetic architecture had negligible effects (Supplementary Figs. 16 and 17 ). Preserving more imaging signals improved MAE and reduced bias. LD regularization controlled the standard error of the estimates (Supplementary Figs. 18 and 19 ). Similar patterns were observed for cross-trait genetic correlation estimates (Supplementary Figs. 20 – 23 ), where we recommend using relatively restrictive LD regularization, such as {75%, 70%} for genotype array data and {90%, 85%} for imputed genotype data, to achieve smaller standard error. We observed high consistency between the mean standard error and empirical standard error of the estimator (Supplementary Fig. 24 ), validating the analytically derived variance. Voxel-level GWAS for brain imaging data in UKB We applied RVGA to analyze hippocampus shape and WM microstructure from 53,454 unrelated European subjects in UKB (Supplementary Data 1 ), and 33,324 of them were used in discovery. We selected the hippocampus shape to show that RVGA can achieve three orders of computational saving, while we utilized the WM microstructure to systematically compare with previous results of IDP-based analysis 9 , 11 on the same data. The shape feature, measured by radial distance from the medial model at each vertex, characterizes morphometric changes perpendicular to the surface. Specifically, the medial model is a mathematical representation of an object’s geometry that is defined by its medial axis or medial surface. WM microstructure is assessed by fractional anisotropy (FA) for each voxel, indicating the restriction of water diffusion along WM tracts. The hippocampus shape is represented as a 3D mesh surface with 30,000 vertices, while the atlas-defined WM tract regions range in size from 88 voxels (inferior fronto-occipital fasciculus) to 3503 voxels (superior longitudinal fasciculus), with a total of 32,217 voxels. To further demonstrate RVGA’s efficiency in whole-cerebral cortex analysis without relying on any atlas for segmentation, we additionally analyzed cortical surface curvature with 59,412 vertices from 28,183 unrelated European subjects ( Supplementary Note ). We adaptively generated LDRs for each WM tract and hippocampus in the left and right hemispheres (each with 15,000 vertices), extracting the top principal components (PCs) that contributed at least 80% of the variance (Methods). The eigenvalues of the hippocampus showed a steeper decay compared to those of WM tracts (Supplementary Fig. 25 ). We constructed 49 LDRs (0.16% of 30,000) to capture 90% of hippocampal signals, whereas 1034 LDRs (3.2% of 32,217) were needed for WM tracts to preserve 80% of the signals (Supplementary Data 2 ). To show that the LDRs captured most of imaging signals, we reconstructed images using the LDRs and bases, calculating a correlation of 0.97 ( s d = 0.03) for the hippocampus and 0.88 ( s d = 0.04) for WM tracts with the raw images. Using imputed genotype data of 7.8 million SNPs (MAF > 0.01), we conducted GWAS on the LDRs. The shared summary statistics for in total 1083 LDRs were 55 GB after processing (Methods), reducing storage burden by 229 times compared to 12,611 GB required to store 62,217 GWAS datasets (the total number of voxels and vertices analyzed) in gzip format. VGWAS summary statistics were reconstructed by saving only significant associations. Specifically, the effective number for each ROI was used to adjust for multiple hypothesis testing (Methods, Supplementary Data 2 ). For the hippocampus (resp. WM tracts), the effective number was 10.1 (resp. 261.8), resulting in a Bonferroni p -value threshold of 5 × 10 −8 /10.1 = 4.94 × 10 −9 (resp. 1.91 × 10 −10 ). We excluded SNPs associated with a small number of voxels by using the wild bootstrap approach 38 (Methods and Supplementary Note ), since the loci derived from these associations are likely due to local perturbations in images. The remaining voxel-variant associations within each ROI were combined into loci using the Peaks algorithm 9 . All significant voxel-variant pairs within a locus were close in genetic distance, with the longest distance from the most significant association being less than 0.25 cM. Compared with combining all voxel-variant associations (815 loci for WM tracts and 134 loci for hippocampus), 37.0% of loci were excluded. For WM tracts, we ended up with 526 loci, with 251 (47.7%) overlapping with significant loci identified in previous studies using the same UKB data and IDP-based approaches 9 , 11 (<0.25 cM), and 49.7% of known loci in these two studies (251 out of 505) were replicated by our study (Fig. 2 A and Supplementary Data 3 ). For the hippocampus, we identified 72 loci, with 33 (45.8%) overlapping with associations in the NHGRI-EBI GWAS catalog and our previous study 13 (Fig. 2 B and Supplementary Data 4 ). Further comparing with results from Hibar et al. 39 where they identified six loci associated with hippocampal volume using data from the ENIGMA Consortium and the CHARGE Consortium ( n = 33,536), four of them were replicated by our study, even if the cohort, image preprocessing pipelines, and features were distinct. In total, we identified 275 novel loci for WM tracts and 39 novel loci for the hippocampus through voxelwise inspection. Fig. 2. The identified genetic loci for WM microstructure and the hippocampus shape. Open in a new tab Each point represents a locus by grouping significant voxel-variant pairs within an ROI such that the longest distance from any variant to the most significant variant was less than 0.25 cM. Each voxel-variant pair can be included in one and only one locus. A previous locus was replicated by our study if the most significant variant in that locus was within 0.25 cM from any of the most significant variants in our loci. We have harmonized the definition of loci in previous studies and the current study before comparison. A Ideogram of loci influencing WM tracts with the p -value threshold 1.91 × 10 −10 (Wald two-sided test, unadjusted p -values), including 275 novel loci and 251 previously identified loci 9 , 11 . B Ideogram of loci influencing hippocampus with the p -value threshold 4.94 × 10 −9 (Wald two-sided test, unadjusted p -values), including 39 novel loci and 33 previously identified loci from our previous study 13 and all hippocampus-related studies on NHGRI-EBI GWAS catalog. We discovered numerous colocalizations with other brain structural and functional measurements, as well as with brain-related phenotypes. At 3q24, the index SNP (i.e., the most significant SNP in a locus), rs2279829, was associated with external capsule (novel) and superior longitudinal fasciculus (previously known) 11 (Fig. 3 A, B). This locus was also linked to insomnia 40 . A novel locus at 13q31.1 was associated with the left hippocampal tail (Fig. 3 E, F) and showed colocalization with smoking 41 – 43 , educational attainment 40 , 44 , and major depressive disorder 45 . Both these two loci can be stringently replicated using an independent dataset (see the next section). We observed that the z-score maps of associations between the index SNP and all voxels/vertices from RVGA were nearly identical to those from VGWAS (Fig. 3 C vs. D, G vs. H). Fig. 3. Selected genetic loci that were associated with subregions in WM tracts/hippocampus, brain structural and functional traits, and brain-related phenotypes. Open in a new tab A At 3q24, we identified a locus with index variant rs2279829 that was simultaneously associated with external capsule (EC, novel, P = 1.79 × 10 −12 ) and superior longitudinal fasciculus (SLF, previously known, P = 4.11 × 10 −16 ). B Significant subregions associated with rs2279829 ( P < 1.91 × 10 −10 , Wald two-sided test, unadjusted p -values) are highlighted in red. C , D z-score maps of associations between the index SNP rs2279829 and all voxels from RVGA and VGWAS, respectively. E At 13q31.1, we discovered a locus with index variant rs76842519 linked to vertices in left hippocampal tail ( P = 5.58 × 10 −12 ). This locus was novel to the hippocampus, while it displayed connections with brain structural and connectivity measurements as well as smoking initiation, educational attainment, and major depressive disorder. F Significant subregions in the left hippocampal tail associated with rs76842519 ( P < 4.94 × 10 −9 , Wald two-sided test, unadjusted p -values) are highlighted in red. G , H z-score maps of associations between the index SNP rs76842519 (left hippocampus) and rs36188842 (right hippocampus, the most significant SNP in the locus) and all vertices from RVGA and VGWAS, respectively. Validating RVGA GWAS results We validated the RVGA results through several complementary analyses. First, we wanted to show that RVGA is aligned well with VGWAS for significant associations. We performed association analyses on 10 randomly selected variants that showed significant associations ( P < 5 × 10 −8 ) in RVGA, along with 90 additional randomly selected variants. We evaluated these variants across all vertices in the left hippocampus and all voxels in the superior fronto-occipital fasciculus using both RVGA and VGWAS (Supplementary Figs. 26 and 27 ). For each variant, we evaluated correlation of effect size/z-score between RVGA and VGWAS. Across all 100 variants, the mean correlation coefficient of genetic effect was 0.90 (sd = 0.05) and the mean correlation coefficient of z-score was 0.84 (sd = 0.07) for the hippocampus, while the mean correlation coefficient of genetic effect was 0.84 (sd = 0.06) and the mean correlation coefficient of z-score was 0.84 (sd = 0.07) for the superior fronto-occipital fasciculus. If focusing on the first 10 variants with significant associations, the counterparts were improved to 0.98 (sd = 0.01) and 0.95 (sd = 0.01) for the hippocampus, and improved to 0.95 (sd = 0.04) and 0.95 (sd = 0.02) for the superior fronto-occipital fasciculus. Moreover, more significant associations were identified by RVGA relative to VGWAS. RVGA aligned well with VGWAS for significant associations, while the noise removal caused the relatively weak alignment for insignificant associations (Discussion). Additionally, we wanted to evaluate RVGA summary statistics across the whole genome. We performed GWAS for the 460,000 genotyped SNPs on 100 randomly selected points from the above two ROIs, comparing RVGA z-scores to those from VGWAS (Supplementary Fig. 28 ). Focusing on SNPs with p -values less than 0.05 in both RVGA and VGWAS, we observed a mean correlation coefficient of 0.99 (sd = 0.005) and a mean RMSE of 0.34 (sd = 0.09) for the left hippocampus. For the superior fronto-occipital fasciculus, the mean correlation and mean RMSE were 0.98 (sd = 0.002) and 0.44 (sd = 0.03), respectively. Considering all SNPs, we observed similar genomic inflation factors ( λ GC ) of 1.05 (sd = 0.01) for RVGA and 1.04 (sd = 0.02) for VGWAS for the left hippocampus, and 1.05 (sd = 0.01) for RVGA and 1.04 (sd = 0.01) for VGWAS for the superior fronto-occipital fasciculus. We further replicated our findings using an independent dataset of 20,130 unrelated European subjects from UKB (Methods and Supplementary Note ). Specifically, we separately constructed FPCA basis for the independent dataset to best capture the image structure. We considered all SNPs in the significant loci from the discovery phase. Because of strong LD among SNPs in a locus, any SNPs being significant in the replication study indicated the locus was replicated. For WM microstructure (Supplementary Data 5 ), we first evaluated alignment of genetic effect size by matching significant voxel-variant pairs in the discovery ( P < 1.91 × 10 −10 ) and in the replication ( P < 0.05/526). The correlation of effect size estimates was 0.96, with 99.3% of effects showing a consistent direction. Then we evaluated loci-level replication. Out of the 526 loci, 513 (97.5%) were replicated (Supplementary Data 3 ), and 369 (70.2%) were replicated additionally considering the effective number ( P < 0.05/526/261.8 = 3.63 × 10 −7 ). For the hippocampus (Supplementary Data 6 ), the correlation of effect size estimates was 0.98, with 99.8% of effects showing a consistent direction. All of the 72 loci were replicated ( P < 0.05/72, Supplementary Data 4 ), and 63 (87.5%) were replicated at the more stringent threshold ( P < 0.05/72/10.1 = 6.88 × 10 −5 ). Focusing on the novel loci, 161 (58.5%) were replicated for WM microstructure and 31 (79.5%) for the hippocampus. These numbers were lower than the overall replication rates because novel loci are associated with local regions, which are generally more difficult to identify and more sensitive to local perturbations such as segmentation or registration errors. Extensive sensitivity analyses, including comparing PCA and FPCA, varying numbers of LDRs, adjusting for additional covariates, and comparing reproducibility between RVGA and VGWAS, are detailed in the Supplementary Note . Heritability and genetic correlation in brain images We estimated voxel heritability and genetic correlations for each ROI using the shared triplets and the LD matrix of 1,160,000 HapMap3 SNPs with regularization of {98%, 95%} (Methods). All ROIs exhibited moderate heritability (Supplementary Data 7 ). The mean heritability for the hippocampus was 18.0% (se = 1.9%), similar to previous results on hippocampal volume (18.8%, se = 1.6%). WM tracts showed a mean heritability of 23.2% (se = 1.9%), ranging from 12.7% for the corticospinal tract to 30.2% for the inferior fronto-occipital fasciculus. We extracted heritability estimates of dMRI TBSS FA traits from Smith et al. 9 and compared them with our mean heritability estimates across voxels for each tract. We took average of heritability estimates for some tracts that were estimated separately for left/right hemispheres in Smith et al.’s paper. The correlation of heritability estimates was 0.56 (two-side t test P = 0.01), and the RMSE was 0.05. Voxelwise heritability estimates revealed significant variations across subregions. Heritability was high (23%) in the CA3 and presubiculum of the hippocampus, while the hippocampal tail, CA1, and subiculum exhibited lower heritability (10%) (Fig. 4 A). In WM tracts, regions of the retrolenticular part of the internal capsule and the superior corona radiata showed high heritability (37%), which decreased to 20% in nearby regions (Fig. 4 A). Fig. 4. Heritability and genetic correlation estimates in brain images. Open in a new tab A The voxelwise heritability of the hippocampus shape and WM microstructure. Refer to Supplementary Data 1 for the full name of WM tracts. B The genetic correlations between the most heritable vertices and other vertices within the left and right hippocampus, respectively. The most heritable vertex is at the subiculum for both the left and right hippocampus. C A comparison of RVGA estimates to SumHer estimates on 100 randomly selected vertices (voxels) on the left hippocampus (the superior fronto-occipital fasciculus). Left column: RVGA implemented FPCA to construct LDRs for a certain proportion of variance, then estimated vertex (voxel) heritability by using LDR summary statistics. Right column: RVGA directly estimated heritability by using summary statistics of each vertex (voxel). D A comparison of genetic correlation estimates between RVGA and LDSC. In each plot, “ r ” refers to the Pearson correlation coefficient, and “ci %” refers to relative confidence interval width, defined as the mean standard error of RVGA estimates to that of SumHer (LDSC) estimates. SumHer (LDSC) was conducted using the summary statistics estimated from VGWAS. We used the “BLD-LDAK” model in SumHer and LD scores estimated by using 1000 Genome data in LDSC. Genetic correlation estimates indicate shared genetic architecture between subregions (Supplementary Data 8 ). Nearby voxels had genetic correlations close to one, but this was not always true for distant voxels. For example, the subiculum of the left hippocampus had a highly positive genetic correlation with CA3 (90%) but a negative correlation with CA1 (−50%). In the right hippocampus, the genetic correlation between the subiculum and CA1 ranged from −45% to 80% (Fig. 4 B). Comparing RVGA heritability estimates with those from SumHer 46 for 100 randomly selected points, we found a correlation coefficient of 0.87 for the left hippocampus and 0.71 for the superior fronto-occipital fasciculus (Fig. 4 C, left column). The relative confidence interval (CI) width between RVGA and SumHer was 0.63 and 0.64, respectively. The two methods were highly consistent when estimating heritability using RVGA directly on each point (Fig. 4 C, right column). Therefore, RVGA did not overestimate heritability but rather enhanced genetic influence by extracting underlying imaging signals through FPCA ( Supplementary Note ). The RVGA genetic correlation estimates closely matched those from LDSC, with correlation coefficients of 0.91 and 0.87, and relative CI widths of 0.56 and 0.64, respectively. RVGA achieved statistically efficient estimates by reducing noise through FPCA and utilizing complete-sample-overlap information. Extensive sensitivity analyses are detailed in the Supplementary Note . Cross-trait genetic correlation between brain images and brain-related phenotypes We gathered summary statistics for 11 complex brain disorders and three brain-related phenotypes (Supplementary Data 9 ), generating atlases of genetic correlations between ROIs and phenotypes (Methods). We applied regularization of {90%, 85%} on the LD matrix in the analysis. To summarize the overall significance of the genetic correlations between ROIs and phenotypes, we used the Cauchy combination strategy 47 to meta-analyze the p -values of all points in an ROI. We controlled the false discovery rate separately at 0.05 for the hippocampus (28 ROI-phenotype pairs, Supplementary Data 10 ) and WM tracts (294 ROI-phenotype pairs, Supplementary Data 11 ), yielding no significant results for the hippocampus and 51 significant results for WM tracts (Fig. 5 A). Of these, 19 (37.3%) corroborated findings from a previous study 11 . Specifically, we replicated widespread genetic correlations for FA of WM tracts with educational attainment, cognitive performance, and intelligence. Additionally, we replicated significant pairs involving the anterior limb of the internal capsule and fornix/stria terminalis with depression, and the superior longitudinal fasciculus with schizophrenia. Comparing with results from Hibar et al. 39 , the authors found a negative genetic correlation (−15.5%, se = 5.3%, P = 0.0034) between hippocampal volume and Alzheimer’s disease. In our results, the hippocampus shape in radial distance also displayed a negative mean genetic correlation across vertices (−3.2%, se = 7.6%), but it was not significant. Fig. 5. Highlighted discoveries from cross-trait genetic correlation analysis between brain images and brain-related phenotypes. Open in a new tab A The overall strength of genetic correlation between ROIs and brain-related phenotypes, which is measured through a p -value generated by meta-analyzing voxelwise genetic correlation p -values ( χ 2 two-sided test) using the Cauchy combination strategy 47 . The asterisks denote significant results by controlling the false discovery rate at 0.05 considering all ROI-phenotype pairs. We conducted multiple hypothesis correction separately for the hippocampus and WM tracts. The boxes denote replicates for the previous study 11 . Italic and underscored phenotypes indicate no evidence of sample overlap with the UKB cohort, while other phenotypes suggest potential sample overlap. Refer to Supplementary Data 1 for the full name of WM tracts. B Selected cross-trait genetic correlations between the ROIs and brain-related phenotypes. EA educational attainment. C , D A comparison of RVGA cross-trait genetic correlation estimates and LDSC counterparts on 100 randomly selected vertices (voxels) on the left hippocampus (the superior fronto-occipital fasciculus). LDSC was conducted using the summary statistics estimated from VGWAS and LD scores estimated by using 1000 Genomes data. In each plot, “ r ” refers to Pearson correlation coefficient, and “ci %” refers to relative confidence interval width, defined as the mean standard error of RVGA estimates to that of LDSC estimates. Inspecting genetic correlations between subregions and phenotypes (Fig. 5 B), we found, for example, that the hippocampal tail, fimbria, CA3, and CA1 had positive correlations (6–10%, se = 3.8%) with educational attainment, while the presubiculum exhibited a negative correlation (−13%, two-sided z test, P = 0.001 < 0.05/10.1). Inconsistent patterns were observed between hemispheres, such as the right anterior corona radiata showing a negative correlation (−9%, P = 4.3 × 10 −5 < 0.05/261.8) with schizophrenia. Moreover, the genu of the corpus callosum was negatively correlated with bipolar disorder, with varying levels across subregions (−16% ~0, se = 2.6%). We validated our findings by comparing RVGA genetic correlation estimates with those from LDSC 48 for 100 randomly selected points. For the left hippocampus and Alzheimer’s disease/educational attainment, the correlation coefficients were 0.78 and 0.73, with relative CI widths of 0.68 and 0.89, respectively (Fig. 5 C). For the superior fronto-occipital fasciculus and the same traits, the coefficients were 0.75 and 0.79, and relative CI widths were 0.63 and 0.93, respectively (Fig. 5 D). Extensive sensitivity analyses are detailed in the Supplementary Note . Computational cost To assess computational efficiency, we benchmarked RVGA using data from the left hippocampus (15,000 vertices and 32,021 subjects) and cortical surface curvature (59,412 vertices and 15,752 subjects), along with 7.8 million variants, on computational clusters equipped with four 2.30 GHz CPUs running in parallel (Table 1 ). For the hippocampus, most steps took less than 30 minutes and used under 25 GB of memory. The most time consuming steps were LDR GWAS and wild bootstrap, taking 4.7 and 5.5 hours, respectively. RVGA demonstrated high efficiency, even for high-resolution images of the whole cerebral cortex. For example, reconstructing whole-genome summary statistics for 59,412 vertices took 32.3 hours, while estimating vertex heritability and genetic correlations between all vertex pairs required only 25 min. Table 1. Benchmarking the computational efficiency for RVGA Time (hours) Memory (GB) Output (GB) FPCA a 1.7 9.7 0.2 Constructing 25 LDRs 0.02 12.0 0.04 LDR GWAS 4.7 20.0 17.0 Processing LDR summary statistics b 0.04 9.5 2.2 Voxel-level summary statistics recovery c 1.2 22.1 0.05 Heritability and genetic correlation 0.03 5.5 1.7 Cross-trait genetic correlation 0.02 3.4 0.002 Wild bootstrap g 5.5 30.0 0.01 FPCA d 2.5 25.5 3.5 Constructing 1750 LDRs 0.3 20.0 0.3 LDR GWAS e 40.8 20.0 411.2 Processing LDR summary statistics b 0.6 20.0 87.2 Voxel-level summary statistics recovery f 32.3 50.0 0.06 Heritability and genetic correlation 0.4 20.0 27.0 Cross-trait genetic correlation 0.3 20.0 0.005 Wild bootstrap g 25.2 60.0 0.01 LD matrix construction h 6.5 4.2 2.8 Open in a new tab The first dataset includes the left hippocampus shape data with 15,000 vertices and 32,021 subjects. The second dataset includes the cortical surface curvature data with 59.412 vertices and 15,752 subjects. A total of 7.8 million common SNPs and 4 CPUs with 2.30 GHz were used in the experiment. a Estimating top 3000 components. b Removing duplicate SNPs, strand-ambiguous SNPs, SNPs with small sample size, missing or infinite values, and indels, resulting in 6.6 million SNPs. Saving the data in HDF5 format of data type float32. c Saving significant associations ( P < 4.94 × 10 −9 ). d Estimating all 15,752 components. e Conducting GWAS for 200 LDRs in a batch. f Saving significant associations ( P < 3.15 × 10 −11 ). g Using 150,000 independent SNPs and computing 50 bootstrap samples. h Using 42,000 subjects with 1,160,000 HapMap3 SNPs from UKB imputed genotype data. Discussion We introduced the novel RVGA framework. RVGA’s major contributions include the ability to demonstrate genetic effects on the human brain at a substantially higher resolution than IDP-based methods, with straightforward potential for extension to other organs using biomedical images. This enables users to precisely identify regions significantly associated with genetic variants, inspect heritability and genetic correlation variations, and evaluate genetic connections between organs and complex disorders. RVGA utilizes FPCA to construct LDRs for images, conducts GWAS on LDRs, and employs functional bases to reconstruct voxel-level summary statistics. This approach significantly reduces computational time, cost, and storage space by orders of magnitude compared to performing GWAS on each individual voxel. Additionally, RVGA optimizes data storage by saving only the summary statistics of LDRs, the functional bases, and the variance-covariance matrix of LDRs. This standard enables the sharing of highly informative voxel-level summary statistics with the community for secondary analyses. In real data analysis, RVGA identified 598 genetic loci, 585 (97.8%) of which can be replicated, and 432 (72.2%) can be replicated at more stringent thresholds. Notably, we discovered novel local genetic connections, such as CA1 and presubiculum with educational attainment, subregions in the anterior corona radiata with schizophrenia, and the genu of the corpus callosum with bipolar disorder. More importantly, we reduced data size for voxel-level secondary analyses by 229 times. We demonstrated RVGA’s robustness and effectiveness through extensive simulations, real data, and sensitivity analyses. When comparing association z-scores between RVGA and VGWAS, we observed that RVGA showed strong alignment with VGWAS for significant associations, although the alignment was weaker for insignificant ones. That is because RVGA effectively smooths images while slightly shrinking each voxel’s value toward zero. In general, this shrinkage leads to a significant reduction in standard error estimates and a small decrease in the absolute effect sizes. For associations with large effect sizes, this reduction can yield smaller p -values. For insignificant associations, standard error estimates are still reduced, but changes in absolute effect sizes are random, as these values are near zero, resulting in weaker alignment. It is the reason why RVGA can boost statistical power while protecting type I error rate. Although top LDRs can capture majority of imaging signals, it is not necessary that no genetic variability is lost. Nevertheless, we have shown that RVGA is more powerful than VGWAS to detect true associations, indicating more noise is removed relative to genetic signals. Additionally, the voxel heritability estimator based on the top LDRs is unbiased under strong noise. Overall, selecting an appropriate number of LDRs can enhance statistical power and computational efficiency. We recommend preserving 80%-90% of imaging signals and ensuring the correlation between the raw and reconstructed images is 0.85-0.95 to balance bias, variance, and computational cost. There are several limitations and potential future directions for RVGA. First, the primary source of bias in RVGA is model mis-specification—it is not able to correctly specify and capture all covariate/genetic effects and imaging signals. Second, RVGA currently restricts GWAS to unrelated subjects and common variants. Future work should address relatedness among subjects 36 and incorporate approaches for whole-genome/exome sequencing data 49 . Third, RVGA is computationally intractable for images with millions of voxels due to the increasing memory demand for FPCA, necessitating more efficient representation learning approaches. Nevertheless, RVGA advances imaging genetic studies to the next level. Beyond traditional biomedical images, RVGA is well-suited to assess other data types, such as microscopy of 3D brains 50 , and transcriptomic and spatial cell-type atlases 51 . Methods Ethics approval The UKB has obtained ethics approval from the Northwest Multi-Centre Research Ethics Committee (MREC, approval number: 11/NW/0382), and obtained written informed consent from all participants prior to the study. Math notations Table 2 demonstrates key math notations used in the article. Table 2. Math notations used in the article Notation Meaning Notation Meaning Y (⋅) Raw image N Number of voxels/vertices X (⋅) Underlying image d Number of variants α (⋅) Fixed covariate effect n Sample size β (⋅) Fixed genetic effect r Number of LDRs η (⋅) Image without covariate and genetic effects ϕ , Φ Functional base ϵ (⋅) White noise λ Image eigenvalue W Covariate matrix u , v Voxel/vertex index Z Genotype matrix ξ , Ξ Low-dimensional representation R LD matrix H kernel bandwidth Ω LD matrix inverse Q (⋅,⋅) Genetic variance/covariance h (⋅) Heritability GC(⋅,⋅) Genetic correlation b k j Projection of β k (⋅) on ϕ j (⋅) ( j -th base) σ 2 Variance of white noise θ i j Projection of η i (⋅) on ϕ j (⋅) ( j -th base) a Polynomial decay rate of image eigenvalues γ Strength of genetic effect size affected by LD α Strength of genetic effect size affected by MAF Open in a new tab A varying coefficient model for imaging genetics In the RVGA framework, the imaging data is assumed to be indexed on a common compact set T ⊂ R 3 for all subjects, capable of capturing curves (1D), surfaces (2D), and volumes (3D) in R 3 . All subjects have N common voxels in the image after appropriate preprocessing. A varying coefficient model for subject i ∈{1,…, n } at voxel v is expressed as: Y i ( v ) = W i ′ ⋅ α ( v ) + Z i ⋅ ′ β ( v ) + η i ( v ) + ϵ i ( v ) . 1 Here, Y i (⋅) represents an N × 1 vector of the image, W i ⋅ is a p × 1 vector of covariates, including the intercept, and α (⋅) is a p × N matrix of fixed coefficients. Z i ⋅ is a d × 1 vector of the genetic profile, and β (⋅) is a d × N matrix of fixed genetic effects. The unexplained image signals are captured by η i (⋅) and ϵ i (⋅). The term η i (⋅) is an N × 1 vector capturing a relatively smooth noise component, while the term ϵ i (⋅), also an N × 1 vector, captures purely random white noise. The genotype matrix Z across all subjects with dimension n × d is assumed to be normalized using the formula Z i k = ( Z i k * − 2 p k ) / 2 p k ( 1 − p k ) . Here, Z i k * is the number of copies of the reference allele for the i -th subject and k -th SNP, and p k is the frequency of the reference allele in the population. The genetic profile Z i ⋅ follows a multivariate normal distribution N ( 0 , R ) , where R = E ( Z i ⋅ Z i ′ ) is a positive-definite LD matrix of d SNPs. Both α l (⋅), l ∈{1,…, p }, and β k (⋅), k ∈{1,…, d }, are fixed smooth functions. The residual term η i (⋅) is modeled as a stochastic process, represented by η i (⋅) ~ SP(0, Σ η (⋅,⋅)). The most common case is Gaussian process; however, the methodology of FPCA does not rely on the Gaussian assumption. For any u , v ∈ T , Cov( η i ( u ), η i ( v )) = Σ η ( u , v ), capturing spatial relationship. The white noise ϵ i ( v ) is assumed to be independent and identically distributed (i.i.d.) across all voxels and subjects. The distribution of ϵ i ( v ) can be normal or other skewed distributions. The covariance of raw image is Cov ( Y i ( u ) , Y i ( v ) ) = Cov ( W i ′ α ( u ) , W i ′ α ( v ) ) + β ( u ) ′ R β ( v ) + Σ η ( u , v ) + σ 2 1 { u = v } , where 1{⋅} is the indicator function. Finally, it is assumed that W i ′ α ( v ) , Z i ′ β ( v ) , η i ( v ), and ϵ i ( v ) are independent of each other. Our varying coefficient model can be regarded as a combination of standard voxel-wise linear regression models and random field theory. At each voxel v , the model decomposes the residual term into individual curve variations ( η i ( v )) and white noise ( ϵ i ( v )). The random field theory is applicable after smoothing Y i ( v ) or reducing ϵ i ( v ) to zero. Compared to voxel-wise linear models, varying coefficient models provide a unified framework for modeling different voxels. We have extensively discussed several key features (e.g., spatial smoothness) of varying coefficient models in Huang et al. 25 . We have used the model to fit almost all MRI imaging phenotypes (e.g., cortical and subcortical structures) used in imaging genetics 52 – 58 , some of which are presented in this paper. Representation learning with functional principal component analysis We apply FPCA to decompose the image Y into smooth imaging signals and white noise. We first smooth the images and then compute a set of smooth functional bases to construct LDRs. In contrast, traditional PCA is less effective at removing noise, as it directly computes bases from raw images. We discussed the differences between FPCA and PCA in dimension reduction and noise removal, and compared RVGA to two related representation learning approaches in the Supplementary Note . A local linear estimator 35 with a Gaussian kernel is used to smooth raw images by approximating the smooth functions underlying the raw images Y . Essentially, each voxel is represented by a weighted average of itself and neighboring voxels. The Gaussian kernel with a specified bandwidth (a hyperparameter) determines the weights and the number of voxels included. We initialize the bandwidth as H = N −1/( d i m +4) , where d i m is the dimension of the image (i.e., 2D or 3D), and use generalized cross-validation (GCV) 34 to select the optimal one from 0.5 H , H , 2 H , 3 H , 5 H , and 10 H . Given the weights, the estimation is equivalent to solving a weighted least squares problem and, therefore, has a closed-form solution. The alternatives include local constant and local quadratic estimators. However, the former suffers from bias and poor boundary behavior, while the latter increases computational burden and estimate variance. Let X i ( v ) = W i ′ α ( v ) + Z i ′ β ( v ) + η i ( v ) denote the real image without white noise, where we assume ∫ T E ( X 2 ) < ∞ , and let the covariance function be C ( u , v ) = Cov( X i ( u ), X i ( v )), which is positive semidefinite. The covariance function admits a spectral decomposition in terms of non-negative eigenvalues λ j that C ( u , v ) = ∑ j = 1 ∞ λ j Φ j ( u ) Φ j ( v ) , u , v ∈ T , 2 where λ 1 ⩾ λ 2 ⩾ ⋯ ⩾ 0, and the functions ϕ 1 , ϕ 2 ,… form an orthonormal basis for the space of all square-integrable functions on T . Denote the sample covariance of the smoothed images as C ^ ( u , v ) . The spectral decomposition is applied to estimate eigenvalues λ ^ 1 ≥ λ ^ 2 ≥ ⋯ , and functional basis ϕ ^ = ( Φ ^ 1 , Φ ^ 2 , … ) such that C ^ ( u , v ) = ∑ j = 1 ∞ λ ^ j Φ ^ j ( u ) Φ ^ j ( v ) (we equivalently use singular value decomposition (SVD) in practice). In real images with discrete voxels, each base is an N × 1 vector. We project Y i onto the subspace spanned by the first r bases and generate LDRs ξ ^ i j = ⟨ Y i , Φ ^ j ⟩ , where j ∈{1, …, r }, r = arg min k ∑ j = 1 k λ ^ j ∑ j = 1 ∞ λ ^ j ≥ q and q is a proportion of variance (⩾80%). The underlying imaging data can be approximately recovered by X i ( ⋅ ) ≈ ∑ j = 1 r ξ ^ i j Φ ^ j ( ⋅ ) . As a guardrail, we suggest the correlation between the raw and reconstructed images greater than 0.85. A higher correlation, ideally 0.90-0.95, is preferable if computationally feasible. For large imaging datasets of high resolution, it may be computationally intractable to compute all λ ’s (PCs) before selecting LDRs. We propose to compute only the top PCs and dynamically select LDRs in such cases (refer to “Subtlety of FPCA” section in Supplementary Note ). Voxel-level GWAS using low-dimensional representations In this section, we first present the association analysis for each voxel-variant pair, then move on to efficient reconstruction using LDR summary statistics. The covariate effects are first removed from the imaging and genetic data using a projection matrix I − M W = I − W ( W ′ W ) − 1 W ′ , where I is the identity matrix, and W is an n × p matrix including the intercept and is assumed to be full rank. Denoting the resulting imaging data as Y ~ = ( I − M W ) Y , the genetic matrix as Z ~ = ( I − M W ) Z . We normalize Z ~ to have a variance of one. A univariate varying coefficient model for VGWAS is Y ~ ( v ) = Z ~ k β k ( v ) + η ( v ) + ϵ ( v ) , 3 where Z ~ k is an n × 1 genotype vector for the k -th SNP, and β k ( v ) is its fixed effect for voxel v . The terms η ( v ) and ϵ ( v ) are n × 1 vectors. Note that Z ~ k ′ Z ~ k = n − 1 after normalization. For each fixed v , model ( 3 ) is essentially a linear model. Therefore, the voxelwise summary statistics can be estimated using the ordinary least square method that β ^ k ( v ) = ( n − 1 ) − 1 Z ~ K ′ Y ~ ( v ) and Var ( β ^ k ( v ) ) = Y ~ ( v ) ′ ( I − Z ~ k Z ~ k ′ n − 1 ) Y ~ ( v ) ( n − 1 ) ( n − p − 1 ) . 4 We use the Wald test to evaluate associations, χ 2 = β ^ k ( v ) 2 Var ( β ^ k ( v ) ) , 5 which follows a χ 1 2 distribution under the null hypothesis that the true marginal effect β k ( v ) = 0, assuming η (⋅) and ϵ (⋅) are Gaussian. When the normal assumption is violated, the test statistic asymptotically follows a χ 1 2 distribution as long as the sample size is large and the variant is common. We use LDR summary statistics to reconstruct β ^ k ( v ) and Var ( β ^ k ( v ) ) . Denote b ^ k ′ = ( n − 1 ) − 1 Z ~ k ′ Ξ ~ r as the genetic effects of SNP k for all r LDRs, and Var ( b ^ k ) = Ξ ~ r ′ ( I − Z ~ k Z ~ r ′ n − 1 ) Ξ ~ r / { ( n − 1 ) ( n − p − 1 ) } as the estimated variance-covariance matrix of LDR genetic effects, where Ξ ~ r = ( I − M W ) Ξ ^ r . Recall that the real imaging data can be approximated by X ( v ) ≈ Ξ ~ r ϕ ^ ( v ) , where ϕ ^ ( v ) is the v -th column of ϕ ^ r ′ . We have β ^ k ( v ) ≈ b ^ k ′ ϕ ^ ( v ) and Var ( β ^ k ( v ) ) ≈ ϕ ^ ( v ) ′ Var ( b ^ k ) ϕ ^ ( v ) . 6 In practice, the genotype matrix is not normalized prior to GWAS, and the actual sample size n k may differ across SNPs due to different missing patterns in genetic profiles. Consequently, the summary statistics and the variance have a more general form that β ^ k ( v ) ≈ ( Z ~ k ′ Z ~ k ) − 1 Z ~ k ′ Ξ ~ r ϕ ^ ( v ) , Var ( β ^ k ( v ) ) ≈ ϕ ^ ( v ) ′ Ξ ~ r ′ Ξ ~ r ϕ ^ ( v ) Z ~ k ′ Z ~ k ( n k − p − 1 ) − ϕ ^ ( v ) ′ b ^ k ′ b ^ k ϕ ^ ( v ) n k − p − 1 , b ^ k j = ( Z ~ k ′ Z ~ k ) − 1 Z ~ k ′ ξ ~ j , Var ( b ^ k j ) = ξ ~ k ′ ξ ~ j Z ~ k ′ Z ~ k ( n k − p − 1 ) − b ^ k j 2 n k − p − 1 . 7 An additional step is to estimate ( Z ~ k ′ Z ~ k ) − 1 from the summary statistics of LDRs that ( Z ~ k ′ Z ~ k ) − 1 = 1 r ∑ j = 1 r ( n k − p − 1 ) Var ( b ^ k j ) + b ^ k j 2 ξ ~ j ′ ξ ~ j . 8 In summary, we only need to compute and store the summary statistics of LDRs, the variance-covariance matrix of the adjusted LDRs 1 n Ξ ~ ′ Ξ ~ r , and the corresponding bases Φ ^ r to recover the VGWAS results. This reduces the storage burden from O ( d N ) to O ( d r + N r + r 2 ), and the computational burden from O ( d N n ) to O ( d N r + d n r ), achieving an average computation reduction by N / r times. Post-GWAS screening based on cluster size The post-GWAS screening is based on evaluating cluster size, defined as the number of associated voxels for a SNP. The key idea is that reliable associations for an SNP should not be restricted to a small number of voxels. We compute the null distribution of cluster size using the wild bootstrap approach 38 . SNPs with a cluster size less than the quantile at 1 − 0.05/(number of loci * effective number) are excluded. The remaining voxel-variant associations are aggregated into loci and reported. Refer to the Supplementary Note for the detailed algorithm of wild bootstrap. Heritability and genetic correlation analysis in images In the heritability and (cross-trait) genetic correlation sections, we ignore the covariate term in model ( 1 ) for the ease of illustration, assuming all covariate effects have been removed. Let the real imaging data without white noise be X i ( v ) = Z i ′ β ( v ) + η i ( v ) , where Z i ⋅ and η i ( v ) are random. The genetic covariance between u , v ∈ T is Q ( u , v ) = Cov ( Z i ′ . β ( u ) , Z i ′ β ( v ) ) = β ( u ) ′ R β ( v ) . 9 The heritability for voxel v is defined as h ( v ) = Q ( v , v ) Var ( X i ( v ) ) = β ( v ) ′ R β ( v ) β ( v ) ′ R β ( v ) + Σ η ( v , v ) , 10 and the genetic correlation between a pair of voxels u , v is defined as G C ( u , v ) = Q ( u , v ) Q ( u , u ) Q ( v , v ) = β ( u ) ′ R β ( v ) β ( u ) ′ R β ( u ) β ( v ) ′ R β ( v ) . 11 Note the variance of white noise σ 2 is not incorporated in the definition of voxel heritability in ( 10 ) because it is introduced by scanner instability rather than by environmental factors. If including σ 2 , the same group of subjects scanned by two image scanners with different levels of white noise would produce inconsistent heritability. Although we cannot observe the real imaging data, the FPCA procedure can reduce the variance of white noise to O ( r / N ) ≈ O (1/ N ), which is negligible compared with β ( v ) ′ R β and Σ η ( v , v ) (assuming they are O (1)). If we were to directly perform analysis on each voxel in the raw imaging data, the heritability estimate would always be downward biased (Fig. 4 C, Supplementary Figs. 14 and 15 , and refer to the Supplementary Note for details). Following the work by Wang et al. 59 , the genetic covariance for u , v ∈ T is estimated using the method of moments that Q ^ ( u , v ) = β ^ ( u ) ′ Ω ^ β ^ ( v ) − tr ( R ^ Ω ^ ) n Cov ( X i ( u ) , X i ( v ) ) , 12 where Ω = R −1 , and tr ( ⋅ ) is the trace of matrix. Check the Supplementary Note for derivation. Both Ω ^ and R ^ can be estimated from reference panels with matched population ancestry. It is recommended to use two different panels without sample overlap to estimate R and Ω for accuracy of tr ( R Ω ) . Cov( X i ( u ), X i ( v )) can be estimated using the variance-covariance matrix of LDRs and the functional bases. A plug-in estimator for heritability is then constructed by h ^ ( v ) = n β ^ ( v ) ′ Ω ^ β ^ ( v ) X ( v ) ′ X ( v ) − tr ( R ^ Ω ^ ) n , Var ( h ^ ( v ) ) = 2 d n 2 + 2 h ( v ) n + 2 h ( v ) 1 − h ( v ) n , 13 and the plug-in estimator for genetic correlation is G C ^ ( u , v ) = Q ^ ( u , v ) Q ^ ( u , u ) Q ^ ( v , v ) , Var ( G C ^ ( u , v ) ) = 4 n + d n 2 1 − G C ( u , v ) 2 2 + 1 n + d n 2 1 − G C ( u , v ) 2 1 h ( u ) + 1 h ( v ) − 2 − 2 G C ( u , v ) 2 Σ η ( u , v ) Q ( u , v ) + d 2 n 2 G C ( u , v ) 2 1 h ( u ) − 1 − Σ η ( u , v ) Q ( u , v ) 2 + d 2 n 2 G C ( u , v ) 2 1 h ( v ) − 1 − Σ η ( u , v ) Q ( u , v ) 2 + d n 2 G C ( u , v ) 4 Σ η ( u , v ) 2 Q ( u , v ) 2 + 1 h ( u ) − 1 1 h ( v ) − 1 − d n 2 G C ( u , v ) 2 Σ η ( u , v ) Q ( u , v ) 1 h ( u ) + 1 h ( v ) − 2 , if Q ( u , v ) ≠ 0 , Var ( G C ^ ( u , v ) ) = d n 2 1 h ( u ) + 1 h ( v ) − 1 + 1 h ( u ) − 1 1 h ( v ) − 1 + 1 n 1 h ( u ) + 1 h ( v ) + 2 , if Q ( u , v ) = 0 . 14 The derivation of the variance is in the Supplementary Note . With the estimate and standard error, we have a χ 2 statistic that asymptotically follows χ 1 2 under the null hypothesis that the true value is 0 59 , 60 . These estimators perfectly adapt to our framework as both β ^ ( u ) ′ Ω ^ β ^ ( v ) and X ( u ) ′ X ( v ) can be first estimated at LDR level and then use the bases to project to the voxel level, thus the high efficiency. The estimators do not depend on the normal distribution assumption, nor an infinitesimal model or specific structure on genetic effects 61 . Therefore, the estimators are expected to be robust to model mis-specification 59 , 62 . Our heritability estimator is related to generalized random effects (GRE) model 62 , and we compared them in the Supplementary Note . Cross-trait genetic correlation analysis between images and non-imaging phenotypes Suppose that we have summary statistics of a non-imaging phenotype U from the identical population, which can be continuous or binary. We aim to evaluate the genetic correlation between U and X ( v ), v ∈ T , i.e., the genetic correlation between U and each voxel in X . Let γ ^ denote the summary statistics of U measured on a common set of SNPs as X ( v ). Suppose that there are n 2 unrelated subjects for U , and n 0 among them are overlapped with the subjects in X ( v ). We estimate the genetic covariance and genetic correlation by Q ^ ( v , γ ) = β ^ ( v ) ′ Ω ^ γ ^ − n 0 tr ( R ^ Ω ^ ) n n 2 Cov ( X i ( v ) , U i ) , G C ^ ( v , γ ) = Q ^ ( v , γ ) Q ^ ( v , v ) Q ^ ( γ , γ ) , Q ^ ( γ , γ ) = γ ^ ′ Ω ^ γ ^ − tr ( R ^ Ω ^ ) n 2 Var ( U i ) , 15 where the subject i was recruited in both studies. Note that in practice, it is unlikely to know the overlapping sample size n 0 and Cov( X i ( v ), U i ) from summary statistics. However, if there is no evidence of subject overlap in two studies, then n 0 = 0, and we have the variance of G C ^ ( v , γ ) in an analytical form 59 that Var ( G C ^ ( v , γ ) ) = G C ( v , γ ) 2 2 h ( v ) 2 d n 2 + G C ( v , γ ) 2 2 h ( γ ) 2 d n 2 2 + 1 h ( v ) h ( γ ) d n n 2 + 1 n 1 − G C ( v , γ ) 2 h ( v ) + 1 n 2 1 − G C ( v , γ ) 2 h ( γ ) , 16 where h ( γ ) = Q ( γ , γ ) Var ( U i ) is the heritability of U . If sample overlap exists between two studies, we implement cross-trait LDSC 48 within RVGA and use the intercept to estimate n 0 Cov( X i ( v ), U i ), as LDSC is exactly a special case of our estimator ( Supplementary Note ). Specifically, we implement cross-trait LDSC for each pair of LDR ξ ~ j and U , obtaining an estimate for n 0 Cov ( ξ ~ j , U i ) . We then use the bases to project n 0 Cov ( ξ ~ j , U i ) onto the voxel level. LD scores are automatically estimated when constructing an LD matrix in RVGA (see the section below). The analytical variance of G C ^ ( v , γ ) is too complicated to derive in this case, and the unknown n 0 and Cov( X i ( v ), U i ) make it impossible to evaluate. We resort to an LD block-wise jackknife with one LD block left out each time to empirically estimate the variance, for which we evenly split the whole genome into approximately 200 blocks of adjacent SNPs. Estimation of LD matrix and LD scores from external datasets The LD matrix R is a d × d positive-semidefinite matrix that measures the correlation among SNPs. The whole genome, excluding the sex chromosomes, can be divided into 1, 703 approximately independent LD blocks (based on the genome of white subjects, GRCh37) 63 . Consequently, we only take into account the correlation among SNPs within each block, assuming that there is no inter-block correlation. In the current version of RVGA, we used only common variants with an MAF greater than 0.01 from the genotype array or HapMap3 SNPs, mainly because it is adequate to capture the majority of heritability. We noticed that regularization on the LD matrix is crucial for accurate and robust estimation. We employ eigen-decomposition and only preserve top eigenvalues. The specific proportion of variance to preserve depends on the type of SNPs included in the LD matrix. More details about LD matrix and LD score estimation are included in the Supplementary Note . p -Value threshold of multiple hypothesis testing across voxels To determine the p -value threshold for multiple hypothesis testing across all voxels in GWAS, heritability and genetic correlation analysis, we resort to the effective number of independent voxels 64 , defined as E = ( 1 + ∑ j = 2 ∞ λ j λ 1 ) 2 1 + ∑ j = 2 ∞ λ j 2 λ 1 2 , 17 where λ j are eigenvalues of Cov( X ( ⋅ )) which are computed in FPCA. For GWAS, we use a Bonferroni threshold of genome-wide significance 5 × 10 −8 / E . And for heritability and genetic correlation, the Bonferroni threshold is 0.05/ E . The pipeline of RVGA The procedure of using RVGA to conduct voxel-level GWAS and secondary analysis includes the following steps. Refer to “The pipeline of RVGA” section in the Supplementary Note for more details. Image loading Images can be in NIFTI, CIFTI, FreeSurfer morphometry data, or text file formats. An additional image or text file for coordinate information is required. Non-imaging phenotypes can also be analyzed using RVGA. FPCA Raw images are initially subjected to kernel smoothing. The optimal bandwidth for smoothing is selected adaptively from a candidate list based on the GCV score. The images are then processed with IncrementalPCA to compute functional bases and eigenvalues. RVGA generates a table summarizing the number of LDRs required to preserve various proportions of image variance, offering guidance on selecting an appropriate number of LDRs. LDR construction Raw images, bases, covariates, and the number of LDRs are input. RVGA constructs the designated number of LDRs and computes the variance-covariance matrix of covariate-effect-removed LDRs. RVGA prints a table of the mean correlation between raw and reconstructed images using varying numbers of LDRs. LDR GWAS LDR GWAS can be internally conducted. GWAS summary statistics processing RVGA processes raw LDR summary statistics and saves in a single file for fast voxel-level GWAS reconstruction and secondary analyses. Specifically, RVGA removes SNPs from LDR or non-imaging phenotype GWAS summary statistics if they exhibit any of the following characteristics: (1) a duplicated rsID; (2) an ambiguous strand; (3) an effective sample size less than 0.67 times the 90th percentile of the sample size; (4) multiple alleles; or (5) a missing or infinite z-score. Voxel-level GWAS reconstruction RVGA is flexible in reconstructing voxel-level summary statistics for different scenarios: (1) scanning the whole genome and all voxels, and saving only significant associations that pass a provided threshold; (2) conducting analysis for the whole genome and a subset of voxels; (3) conducting analysis for selected variants or a genome segment across all voxels. Post-GWAS screening RVGA computes an empirical null distribution of cluster size, defined as the number of associated voxels for a SNP, using the wild bootstrap approach (refer to the Supplementary Note for specific algorithm). Users can manually exclude SNPs with a small cluster size based on the quantile at level 1 − 0.05/(number of loci * effective number). The remaining voxel-variant associations can be aggregated into loci and reported. LD matrix estimation RVGA estimates the LD matrix and its inverse from a pair of PLINK bfiles based on a specified regularization level. Heritability and (cross-trait) genetic correlation analysis The inputs include processed LDR summary statistics, processed non-imaging phenotype summary statistics (optional), the LD matrix and its inverse, bases, and the variance-covariance matrix of covariate-effect-removed LDRs. Computational complexity For an imaging-genetic study with N total voxels, r LDRs, p covariates, n subjects, n 1 external subjects, m LD blocks, and d SNPs, the computational complexities are as follows: Adaptively selecting the optimal bandwidth and estimating a sparse smoothing matrix using the local-linear method takes the maximum O ( N 2 + n N ) time and O ( N 2 + n N ) memory. Randomized SVD for the smoothed imaging data takes the maximum O ( n N min { n , N } ) time and O ( N 2 + n N ) memory. Constructing LDRs takes O ( N n r ) time; computing the variance-covariance matrix of LDRs takes O ( r 2 n + r p n + p 3 + r 2 p ) time. Storing the data (including the top r bases) takes O ( n r + N r + r 2 ) storage space. Computing summary statistics for LDRs takes O ( d n r ) time and O ( d r ) storage space. Assume each block in the LD matrix has O ( d / m ) SNPs, then estimating the whole LD matrix takes O ( n 1 d 2 / m ) time and O ( d 2 / m ) storage space. Computing heritability and genetic correlation takes O ( r d 2 / m + N 2 r ) time. Specifically, B ^ ′ Ω ^ B ^ takes O ( r d 2 / m + r 2 d ) time, tr ( R ^ Ω ^ ) takes O ( d 2 / m ) time, and mapping the results on the low-dimensional space to the original space takes O ( N r 2 + N 2 r ) time. Saving heritability and genetic correlation estimates takes O ( N ) and O ( N 2 ) storage space, respectively. Cumulatively, the proposed method takes the maximum O ( n N min { n , N } + d n r + n 1 d 2 / m ) time, O ( N 2 + n N ) memory, and O ( N 2 + d r + d 2 / m ) storage space. In the “The pipeline of RVGA” section in the Supplementary Note , we outline the input and output of each step in the pipeline, along with general strategies to enhance computational efficiency. Imaging acquisition and feature generation We obtained structural MRI (sMRI) and diffusion MRI (dMRI) imaging data through UK Biobank application 22783. The UKB team has already preprocessed and quality-controlled the data before releasing it. For detailed information on image acquisition and preprocessing, refer to the protocol ( https://biobank.ctsu.ox.ac.uk/crystal/crystal/docs/brain_mri.pdf ). Using FMRIB Software Library (FSL v5.0.9) ( https://fsl.fmrib.ox.ac.uk/ ), sMRI were linearly registered into a standard brain space (MNI152). We segmented the left and right hippocampus from the sMRI. Surface meshes were constructed by using the topology-preserving level set method 65 and the marching cube algorithm 66 , parameterized with refined triangular meshes using the topological optimization algorithm 67 and holomorphic flow segmentation method 68 . Then the images were registered to a common rectangular grid template using the surface fluid registration algorithm 69 . Radial distance was computed for each vertex. The radial distance represents the Euclidean distance between each vertex and the medial axis, the geometric center of the isoparametric curve. Both the left and right hippocampus images had a resolution of 100 × 150 × 1 (15,000 vertices). The pipeline is available at https://www.nitrc.org/frs/?group_id=1461 . The detailed steps can be found in the Supplementary Note of ref. 13 . Regarding dMRI features, each individual fractional anisotropy (FA) image was generated by fitting diffusion tensor imaging (DTI) models using the FSL software. Our analysis employed the ENIGMA-DTI pipeline, a standardized set of procedures for processing DTI data ( http://enigma.ini.usc.edu/protocols/dti-protocols/ ). We first linearly registered each FA image to the ENIGMA FA template, which was set at a 1 × 1 × 1mm 3 spatial resolution in the MNI152 standard space. We then applied nonlinear registration techniques to further refine the alignment of the FA images to the standard space and masked the registered FA images using the template brain mask. We manually checked the registration performance and removed those with low quality. The ENIGMA skeleton, representing the major WM pathways, was projected onto the registered images. We finally extracted 21 predefined WM tracts according to the JHU ICBM-DTI-81 WM atlas. Tract resolutions ranged from 88 (inferior fronto-occipital fasciculus) to 3503 voxels (superior longitudinal fasciculus). More details for implementing the pipeline can be found in ref. 11 . Real data analysis We collected the UKB phases 1–3 sMRI and dMRI data for 35,000 subjects of European ancestry. To remove related subjects, we computed genetic relatedness matrix using the genotype array data of autosomes with 460,000 common SNPs (MAF > 0.01) through GCTA 70 (v1.93.2beta, https://yanglab.westlake.edu.cn/software/gcta/#Overview ). We removed one of a pair of related subjects with the threshold 0.05 (–grm-cutoff 0.05), which excluded the fourth and more related relatives. We further excluded subjects with excessive heterozygosity (Field ID 22027), discrepancies between reported and genetic gender (Field ID 22001), potential sex chromosome anomalies (Field ID 22019). The final sample size was 33,324. The LDRs were constructed to capture more than 80% of the variance for each WM tract. For hippocampus, we captured 90% of variance due to rapid decay of eigenvalues. The LDRs were adjusted for the intercept, age, sex, age 2 , age × sex, age 2 × sex, assessment center (Data Field 54), and 40 genetic PCs to estimate the variance-covariance matrix. Imputed genotype data was used in LDR GWAS. During quality control, we removed SNPs with (1) an MAF less than 0.01; (2) a p -value of Hardy–Weinberg equilibrium (HWE) test less than 10 −7 ; (3) a missing call rate greater than 10%; and (4) an imputation score less than 0.9, resulting in approximately 7.8 million autosomal SNPs. After processing LDR GWAS summary statistics, we reconstructed voxel-level summary statistics and saved only significant signals. The significance threshold was determined by the effective number (Methods, Supplementary Data 2 ). For hippocampus, we considered left and right hemispheres altogether, and the significance threshold was 4.94 × 10 −9 . We considered all WM tracts altogether, and the threshold was 1.91 × 10 −10 . We extracted 150,000 independent SNPs with a pairwise LD less than 0.1 from 460,000 genotyped SNPs in UKB and computed a null distribution of cluster size using the wild bootstrap approach with 50 bootstrap samples, resulting 7,500,000 points (Methods and Supplementary Note ). SNPs with a cluster size less than the quantile at 1 − 0.05/(number of loci * effective number) were excluded. Note the number of loci can be determined by aggregating all voxel-variant associations before screening. We applied the Peaks algorithms 9 to group the remaining voxel-level associations within each ROI. Specifically, a locus is defined as a collection of significant voxel-variant pairs such that the distance from any variant in the locus to the most significant variant is less than 0.25 cM. Any voxel-variant pair can be included in one and only one locus. A locus in previous studies was replicated by our study if the most significant variant in that locus was within 0.25 cM from any of the most significant variants in our loci. For WM tracts, we compared with previous studies using UKB WM FA traits 9 , 11 . For the hippocampus, we checked our previous study 13 and all the hippocampus-related associations published on the NHGRI-EBI GWAS catalog (up to April 2024). We estimated LD matrices including genotyped SNPs or imputed HapMap3 SNPs. For genetyped SNPs, two datasets of unrelated white subjects from UKB were extracted, each with a sample size of 8400 and no sample overlap. SNPs in the major histocompatibility complex (MHC, GRCh37: chr6 28.5M–33.5M) and/or those with an MAF less than 0.01 and/or a p -value of HWE test less than 10 −7 and/or a missing call rate greater than 10% were removed, resulting in approximately 460,000 SNPs. The genome was partitioned into 1703 approximately independent LD blocks 63 . Missing values in a SNP vector were imputed by the sample mean. The LD matrix for each block was first estimated by sample correlation, then subject to eigen decomposition, and 85% and 80% of variance were preserved for the LD matrix and its inverse, respectively. To enhance statistical efficiency in cross-trait genetic correlation estimates, we used a more restrictive regularization of {75%, 70%}. For imputed HapMap3 SNPs, we used two datasets each with a sample size of 42, 000. The SNPs were screened using the same procedure, resulting in 1, 160, 000 SNPs in analysis. The regularization was {98%, 95%} for heritability and genetic correlation and {90%, 85%} for cross-trait genetic correlation. We gathered summary statistics for 11 complex brain disorders, including Alzheimer’s disease 71 , schizophrenia 72 , insomnia 40 , autism spectrum disorder (ASD) 73 , bipolar disorder (BD) 74 , major depressive disorder (MDD) 75 , attention-deficit hyperactivity disorder (ADHD) 76 , neuroticism, neuroticism subclusters (depressed affect and worry), and depression 77 . Additionally, we included three brain-related phenotypes: educational attainment 44 , cognitive performance 44 , and intelligence 78 . Detailed information for each dataset is provided in Supplementary Data 9 . Most summary statistics were generated through meta-analysis on multiple cohorts and consortiums. For those involving UKB cohorts, we accounted for potential sample overlap in RVGA. Validation of the RVGA voxel-level GWAS results We compared voxel-level summary statistics produced by RVGA to those from VGWAS (PLINK2) using correlation and RMSE of beta coefficients and z-scores. In the first phase, we randomly selected 100 variants for analysis in the left hippocampus and another 100 variants for the superior fronto-occipital fasciculus, where the first 10 variants in each ROI were randomly selected from significant loci. We used a mixture of variants with or without significant associations to mimic the reality that most of variants are insignificant. Another reason might be minor, but we wanted to visualize that RVGA can identify more significant voxel-variant associations than VGWAS in real data, so we randomly selected 10 variants with significant associations and visualized all of them in one figure (Supplementary Figs. 26 and 27 ). In the second phase, we randomly selected 100 points from the above two ROIs and scanned the whole genome using 460,000 genotyped SNPs. We collected UKB phases 4–6 sMRI and dMRI data, including 20,130 unrelated subjects of European ancestry, for replication. The imaging and genotype data were processed using the same pipelines and thresholds as before. We independently estimated functional bases while adjusted for the identical list of covariates. Specifically, we extracted all SNPs in significant loci from the discovery phase and evaluated them in the replication study. Any SNPs in a locus being significant in the replication study indicated the locus was replicated. We used two thresholds in replication: 0.05/(number of loci) and 0.05/(number of loci * effective number). Refer to the Supplementary Note for a brief review of replication studies in the literature, and more details about the replication study. Evaluation of the RVGA heritability and genetic correlation estimates We compared RVGA heritability estimates to those from SumHer 46 . We downloaded pre-computed taggings of white subjects produced by using UKB genotyped SNPs and HapMap3 SNPs, respectively, corresponding to the “BLD-LDAK” model ( https://dougspeed.com/pre-computed-tagging-files/ ). RVGA genetic correlation and cross-trait genetic correlation estimates were compared to those from LDSC 48 . For genotyped SNPs, we estimated LD scores using a genotype array dataset of sample size 8400 through the LDSC command line tool ( https://github.com/bulik/ldsc?tab=readme-ov-file#ldsc-ld-score-v101 ), setting LD window as 10,000 Kb (–ld-wind-kb 10,000). For HapMap3 SNPs, we downloaded pre-computed LD scores from the 1000 Genomes data. SumHer and LDSC were directly applied to each voxel using the summary statistics from VGWAS. In cross-trait genetic correlation analysis, if the cohort in non-imaging phenotypes had no evidence of sample overlap with the UKB cohort, we employed the RVGA estimator without correction, and constrained the intercept in cross-trait LDSC as 0 (–intercept-gencov 0,0). Otherwise, we employed the RVGA estimator considering sample overlap, and did not constrain the LDSC intercept. We calculated Pearson correlation coefficient and the relative confidence interval width, defined as the mean standard error of RVGA estimates to that of SumHer/LDSC estimates across 100 points. Simulation studies We employed real genotype array data or imputed genotype data (HapMap3 SNPs) from UKB in the simulation. We randomly extracted n = 10,000 unrelated white subjects and all common variants (MAF > 0.01) on chromosome 10 (23,203 genotyped SNPs) or on chromosome 19 (21,576 imputed SNPs). We standardized each SNP by the formula Z i k = ( Z i k * − 2 p k ) / 2 p k ( 1 − p k ) , where Z i k * ∈ { 0 , 1 , 2 } is the number of reference alleles, and p k is the sample MAF. To simulate imaging data, we employed a varying coefficient model Y i ( v ) = W i ′ α ( v ) + Z i ′ . β ( v ) + η i ( v ) + ϵ i ( v ) , X i ( v ) = W i ′ . α ( v ) + Z i ′ . β(v) + η i ( v ) , 18 where i ∈ {1, …, n } and v ∈ { 1 N , 2 N , … , 1 } . Moreover, W i ⋅ is a 2 × 1 covariate vector, where the first covariate was sampled from B e r n o u l l i (0.5) and the second one was sampled from N ( 0 , 0.01 ) . Each fixed covariate effect α k ( v ) = ∑ j = 1 200 α 0 k j Φ j ( v ) , where α 0 k j ~ N ( 0 , ( j + 3 ) − 3 ) for k = 1, 2. A set of SNPs C was randomly chosen as causal SNPs. For k ∈ C , β k ( v ) = ∑ j = 1 200 b k j Φ j ( v ) , with Φ j ( v ) = 2 cos { ( j − 1 ) π v } and b k j ~ N 0 , 1 ∣ C ∣ ( j + 3 ) − 1.2 a l k − γ { 2 p k ( 1 − p k ) } 1 + α , where ∣ C ∣ is the number of causal SNPs; a > 1 is a polynomial decay rate; γ ∈ {1, 0} specifies whether effect sizes are related to the LD weights 1/ l k ; α ∈ {−1, −0.25} indicates strong/weak inverse relationship between MAF and effect sizes. For k ∉ C , β k (⋅) = 0. The underlying imaging data without white noise was decomposed as X i ( v ) = ∑ j = 1 200 ξ i j Φ j ( v ) with ξ i j ~ N ( 0 , λ j ) , where λ j = 2( j +2) − a . The non-genetic effect η i ( v ) = ∑ j = 1 200 θ i j Φ j ( v ) , where θ i j ~ N ( 0 , max { λ j − Var ( W i ′ α 0 j + Z i ′ b j ) , 0 } ) . We adjusted the genetic variance of each voxel to achieve a fixed heritability level. Finally, we generated ϵ i ( v ) following N ( 0 , σ 2 ) (normal noise) or Rayleigh distribution with parameter σ , where σ 2 was adaptively set for a specific noise percentage w = σ 2 /(Var( X (⋅)) + σ 2 ). For assessing type I error rate, we used the reduced model without genetic effect Y i ( v ) = W i ′ α ( v ) + η i ( v ) + ϵ i ( v ) and set N = 200, n = 10,000, a = {1.8, 2.5}, w = {50%, 20%, 0}, and the proportion of variance to keep with LDRs q = {0.75, 0.80, 0.85, 0.90, 0.95}. For each scenario, 1000 replicates were generated. We investigated the following scenarios for statistical power: N = 200, n = 10,000, ∣ C ∣ / d = { 1 % , 20 % } (proportion of causal SNPs or polygenicity), h = {0.03, 0.1, 0.3}, α = −1, γ = 0, a = {1.8, 2.5}, w = {50%, 20%, 0}, noise following N ( 0 , σ 2 ) or Rayleigh distribution with scale σ , and the proportion of variance to keep with LDRs q = {0.75, 0.80, 0.85, 0.90, 0.95}. For each scenario, 100 replicates were generated. For evaluating heritability and genetic correlation estimates, we considered the following scenarios: N = 200, n = 10,000, ∣ C ∣ / d = { 1 % , 20 % } , h = {0.03, 0.1, 0.3} ( h = 0.03 only for heritability as it is too small to get stable genetic correlation estimates), α = {−1, −0.25}, γ = {1, 0}, a = 1.8, w = 20%, normal noise, and the proportion of variance to keep with LDRs p r o p = {0.75, 0.80, 0.85, 0.90, 0.95}. For each scenario, 100 replicates were generated. We further simulated a single trait U to evaluate cross-trait genetic correlation estimates. The true genetic effect β U was simulated by β U = 0.3 β ( N / 2 ) Var ( Z i ′ β ( N / 2 ) ) , where β ( N /2) means the genetic effect for the N /2-th voxel, and the observed data U = Z U β U + ϵ U was generated, where Z U was extracted from the real data with no sample overlap with Z and ϵ U ~ N ( 0 , Var ( Z U , i ′ . β U ) I ) . That indicates the heritability of U was fixed at 0.5. The LD matrix and its inverse were independently generated, each using 10, 000 unrelated white subjects. For genotyped SNPs, four regularization levels for the LD matrix and its inverse were assessed: {0.90, 0.85}, {0.85, 0.80}, {0.80, 0.75}, and {0.75, 0.70}. For imputed SNPs, the setup was {0.98, 0.95}, {0.95, 0.90}, {0.90, 0.85}, and {0.85, 0.80}. The type I error was evaluated by the proportion of null tests with a p -value less than 10 −2 , 10 −3 , and 10 −4 . The statistical power was quantified by the proportion of causal variants identified at level 10 −4 . The distance between RVGA summary statistics and those from VGWAS was evaluated by RMSE. For heritability and (cross-trait) genetic correlation analysis, we employed MAE to compare the estimates with the ground truth. We additionally compared the average estimate to the average true value across all voxels. Reporting summary Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article. Supplementary information Supplementary Information (30.7MB, pdf) 41467_2026_69816_MOESM2_ESM.pdf (118.1KB, pdf) Description of Additional Supplementary Files Supplementary Dataset 1-15 (890.8KB, xlsx) Reporting Summary (4MB, pdf) Transparent Peer Review File (19.3KB, pdf) Acknowledgements We thank research assistants (X. Qian, S. Huang, P. Guan, J. Chen, M. Hu, X. Qi, B. Tang, and V. Pathare) in the UNC BIGS2 team for downloading and preprocessing raw images. We thank the individuals represented in the UKB study for their participation and the research teams for their work in collecting, processing, and disseminating these datasets for analysis. We thank University of North Carolina at Chapel Hill and the Research Computing groups for providing computational resources and support that have contributed to the research results. This research has been conducted using the UK Biobank resource (application number 22783), subject to a data transfer agreement. We thank S. Gao for helping polish the figures. We used imaging icons from BioRender.com in Fig. 1 . This work was partially supported by the National Institute on Aging (NIA) of the National Institutes of Health (NIH) grants [R01AG085581, RF1AG082938, U01AG088667, R01AG098697 to T.L. and H.Z.], NIH grants [U01HG011720 and R01MH125236 to Y.L., 1U01AG088667-01 to J.S., R01MH136055 to H.Z. and T.L.; K01AG095286 and R21HD120911 to T.L.], and National Institute of Child Health and Human Development grant [P50HD103573 to Y.L.]. Additionally, this work was partially supported by the National Science Foundation (NSF) grants [DMS-2346292 and DMS-2434666 to E.F.]. The content is solely the responsibility of the authors and does not necessarily represent the official views of these institutes. Author contributions Z.J. and H.Z. proposed the concept and designed the methodology. H.Z. supervised the project. Z.J. performed the statistical analysis, visualized the results, developed and tested the software. T.L. designed the pipeline for preprocessing raw images. Z.J. drafted the initial manuscript. Z.J., J.S., P.S., E.F., Y.L., and H.Z. suggested revision ideas and revised the manuscript. All authors critically reviewed the manuscript and approved the final version. Peer review Peer review information : Nature Communications thanks the anonymous reviewers for their contribution to the peer review of this work. A peer review file is available. Data availability The triplets of summary statistics, as well as the LD matrix (including the pre-computed LD scores in the last column in the. ldinfo file) generated in this study, have been deposited in Zenodo under accession code 13787684 79 , 10.5281/zenodo.13787684. The raw structural MRI and diffusion MRI imaging data, as well as genetic data is under restricted access for privacy issues; access can be obtained by application at https://www.ukbiobank.ac.uk . The data used in this study were obtained from UK Biobank application 22783. The colocalization results were extracted from the NHGRI-EBI GWAS catalog, https://www.ebi.ac.uk/gwas/ . The pre-computed LD scores for the Hapmap3 SNPs are at https://console.cloud.google.com/storage/browser/broad-alkesgroup-public-requester-pays/LDSCORE?pageState=(%22StorageObjectListTable%22:(%22f%22:%22%255B%255D%22)) . Links to summary statistics: schizophrenia, https://figshare.com/articles/dataset/scz2022/19426775?file=34517828 ; Alzheimer’s disease, https://ctg.cncr.nl/documents/p1651/PGCALZ2sumstatsExcluding23andMe.txt.gz ; insomnia, https://ctg.cncr.nl/documents/p1651/insomnia_ukb2b_EUR_sumstats_20190311_with_chrX_mac_100.txt.gz ; educational attainment, https://thessgac.com/papers/3 ; cognitive performance, https://thessgac.com/papers/3 ; autism spectrum disorder, https://pgc.unc.edu/for-researchers/download-results/ ; bipolar disorder https://pgc.unc.edu/for-researchers/download-results/ ; major depressive disorder, https://pgc.unc.edu/for-researchers/download-results/ ; attention deficit hyperactivity disorder, https://pgc.unc.edu/for-researchers/download-results/ ; intelligence, https://ctg.cncr.nl/documents/p1651/SavageJansen_IntMeta_sumstats.zip ; neuroticism, https://ctg.cncr.nl/documents/p1651/sumstats_neuroticism_ctg_format.txt.gz ; depression, https://ctg.cncr.nl/documents/p1651/sumstats_depression_ctg_format.txt.gz ; depression affect subcluster, https://ctg.cncr.nl/documents/p1651/sumstats_depressed_affect_ctg_format.txt.gz ; worry subcluster, https://ctg.cncr.nl/documents/p1651/sumstats_worry_ctg_format.txt.gz . All other data supporting the findings of this study are available in the article and the Supplementary Information files. Code availability RVGA is a component in the Highly Efficient Imaging Genetics (HEIG) toolbox, which is a comprehensive solution for voxel-level imaging genetic analysis. HEIG 80 (v1.6.1-alpha) is an open-source Python software, https://github.com/Zhiwen-Owen-Jiang/HEIG (citable version https://zenodo.org/records/18049388 ) with tutorial, https://github.com/Zhiwen-Owen-Jiang/HEIG/wiki and example data, 10.5281/zenodo.14214075; The example code for reproducing the main results, https://github.com/Zhiwen-Owen-Jiang/HEIG/tree/pub/misc/code_manuscript/code.md ; LDSC (v1.0.1), https://github.com/bulik/ldsc ; SumHer, https://dougspeed.com/sumher/ ; GCTA (v1.93.2beta), https://yanglab.westlake.edu.cn/software/gcta/#Overview ; PLINK2 (v2.00a3LM), https://www.cog-genomics.org/plink/2.0/ ); Peaks algorithm (v1.0), https://github.com/wnfldchen/peaks . Competing interests P.S. reports the following potentially competing financial interests (2024): Neumora Therapeutics (advisory committee, shareholder). To the best of his knowledge, these are unrelated to this paper/project. The remaining authors declare no competing interests. Footnotes Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. Supplementary information The online version contains supplementary material available at 10.1038/s41467-026-69816-z. References 1. Scangos, K. W., State, M. W., Miller, A. H., Baker, J. T. & Williams, L. M. New and emerging approaches to treat psychiatric disorders. Nat. Med. 29 , 317–333 (2023). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 2. Moreau, C. A. et al. Dissecting autism and schizophrenia through neuroimaging genomics. Brain 144 , 1943–1957 (2021). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 3. Bearden, C. E. & Thompson, P. M. Emerging global initiatives in neurogenetics: the enhancing neuroimaging genetics through meta-analysis (ENIGMA) consortium. Neuron 94 , 232–236 (2017). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 4. Grasby, K. L. et al. The genetic architecture of the human cerebral cortex. Science 367 , eaay6690 (2020). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 5. Jack, C. R. et al. Tracking pathophysiological processes in alzheimer’s disease: an updated hypothetical model of dynamic biomarkers. Lancet Neurol. 12 , 207–216 (2013). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 6. Jack, C. R. et al. Hypothetical model of dynamic biomarkers of the alzheimer’s pathological cascade. Lancet Neurol. 9 , 119–128 (2010). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 7. Jiang, Z. et al. The X chromosome’s influences on the human brain. Sci. Adv. 11 , eadq5360 (2025). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 8. Elliott, L. T. et al. Genome-wide association studies of brain imaging phenotypes in UK Biobank. Nature 562 , 210–216 (2018). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 9. Smith, S. M. et al. An expanded set of genome-wide association studies of brain imaging phenotypes in UK Biobank. Nat. Neurosci. 24 , 737–745 (2021). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 10. Zhao, B. et al. Genome-wide association analysis of 19,629 individuals identifies variants influencing regional brain volumes and refines their genetic co-architecture with cognitive and mental health traits. Nat. Genet. 51 , 1637–1644 (2019). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 11. Zhao, B. et al. Common genetic variation influencing human white matter microstructure. Science 372 , eabf3736 (2021). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 12. Zhao, B. et al. Common variants contribute to intrinsic human brain functional networks. Nat. Genet. 54 , 508–517 (2022). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 13. Zhao, B. et al. Genetic influences on the shape of brain ventricular and subcortical structures. Preprint at medRxiv 2022–09 10.1101/2022.09.26.22279691 (2022). 14. Stein, J. L. et al. Identification of common variants associated with human hippocampal and intracranial volumes. Nat. Genet. 44 , 552–561 (2012). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 15. Patel, K. et al. Unsupervised deep representation learning enables phenotype discovery for genetic association studies of brain imaging. Commun. Biol. 7 , 414 (2024). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 16. Wen, J. et al. Genomic loci influence patterns of structural covariance in the human brain. Proc. Natl. Acad. Sci. USA 120 , e2300842120 (2023). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 17. Zhao, B. et al. Eye-brain connections revealed by multimodal retinal and brain imaging genetics. Nat. Commun. 15 , 6064 (2024). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 18. Stein, J. L. et al. Voxelwise genome-wide association study (VGWAS). Neuroimage 53 , 1160–1174 (2010). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 19. Hibar, D. P. et al. Voxelwise gene-wide association study (vgenewas): multivariate gene-based association testing in 731 elderly subjects. Neuroimage 56 , 1875–1891 (2011). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 20. Nichols, T. & Hayasaka, S. Controlling the familywise error rate in functional neuroimaging: a comparative review. Stat. Methods Med. Res. 12 , 419–446 (2003). [ DOI ] [ PubMed ] [ Google Scholar ] 21. Silver, M. et al. False positives in neuroimaging genetics using voxel-based morphometry data. Neuroimage 54 , 992–1000 (2011). [ DOI ] [ PMC free article ] [ PubMed ] 22. Hua, W.-Y., Nichols, T. E., Ghosh, D. & Initiative, A. D. N. Multiple comparison procedures for neuroimaging genomewide association studies. Biostatistics 16 , 17–30 (2015). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 23. Ge, T., Feng, J., Hibar, D. P., Thompson, P. M. & Nichols, T. E. Increasing power for voxel-wise genome-wide association studies: the random field theory, least square kernel machines and fast permutation procedures. Neuroimage 63 , 858–873 (2012). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 24. Huang, M. et al. Fvgwas: Fast voxelwise genome wide association analysis of large-scale imaging genetic data. Neuroimage 118 , 613–627 (2015). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 25. Huang, C. et al. Fgwas: functional genome wide association analysis. Neuroimage 159 , 107–121 (2017). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 26. Wen, C., Ba, H., Pan, W. & Huang, M. Co-sparse reduced-rank regression for association analysis between imaging phenotypes and genetic variants. Bioinformatics 36 , 5214–5222 (2020). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 27. Zhu, X., Zhang, W. & Fan, Y. A robust reduced rank graph regression method for neuroimaging genetic analysis. Neuroinformatics 16 , 351–361 (2018). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 28. Huang, M. et al. Spatial correlations exploitation based on nonlocal voxel-wise GWAS for biomarker detection of ad. NeuroImage Clin. 21 , 101642 (2019). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 29. Li, J., Wang, Z., Li, R. & Wu, R. Bayesian group lasso for nonparametric varying-coefficient models with application to functional genome-wide association studies. Ann. Appl. Stat. 9 , 640 (2015). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 30. Lu, Z.-H. et al. Bayesian longitudinal low-rank regression models for imaging genetic data from longitudinal studies. NeuroImage 149 , 305–322 (2017). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 31. Stingo, F. C., Guindani, M., Vannucci, M. & Calhoun, V. D. An integrative Bayesian modeling approach to imaging genetics. J. Am. Stat. Assoc. 108 , 876–891 (2013). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 32. Fan, C. C. et al. Multivariate genome-wide association study on tissue-sensitive diffusion metrics highlights pathways that shape the human brain. Nat. Commun. 13 , 2423 (2022). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 33. Zhou, H., Yao, F. & Zhang, H. Functional linear regression for discretely observed data: from ideal to reality. Biometrika 110 , 381–393 (2023). [ Google Scholar ] 34. Zhu, H., Li, R. & Kong, L. Multivariate varying coefficient model for functional responses. Ann. Stat. 40 , 2634 (2012). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 35. Fan, J. Local Polynomial Modelling and its Applications: Monographs on Statistics and Applied Probability 66 (Routledge, 2018). 36. Mbatchou, J. et al. Computationally efficient whole-genome regression for quantitative and binary traits. Nat. Genet. 53 , 1097–1103 (2021). [ DOI ] [ PubMed ] [ Google Scholar ] 37. Speed, D., Hemani, G., Johnson, M. R. & Balding, D. J. Improved heritability estimation from genome-wide snps. Am. J. Hum. Genet. 91 , 1011–1021 (2012). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 38. Zhu, H. et al. Frats: Functional regression analysis of dti tract statistics. IEEE Trans. Med. Imaging 29 , 1039–1049 (2010). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 39. Hibar, D. P. et al. Novel genetic loci associated with hippocampal volume. Nat. Commun. 8 , 13624 (2017). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 40. Watanabe, K. et al. Genome-wide meta-analysis of insomnia prioritizes genes associated with metabolic and psychiatric pathways. Nat. Genet. 54 , 1125–1132 (2022). [ DOI ] [ PubMed ] [ Google Scholar ] 41. Saunders, G. R. et al. Genetic diversity fuels gene discovery for tobacco and alcohol use. Nature 612 , 720–724 (2022). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 42. Liu, M. et al. Association studies of up to 1.2 million individuals yield new insights into the genetic etiology of tobacco and alcohol use. Nat. Genet. 51 , 237–244 (2019). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 43. Pasman, J. A. et al. Genetic risk for smoking: disentangling interplay between genes and socioeconomic status. Behav. Genet. 52 , 92–107 (2022). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 44. Lee, J. J. et al. Gene discovery and polygenic prediction from a genome-wide association study of educational attainment in 1.1 million individuals. Nat. Genet. 50 , 1112–1121 (2018). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 45. Investigators, G., Investigators, M. & Investigators, S. D. Common genetic variation and antidepressant efficacy in major depressive disorder: a meta-analysis of three genome-wide pharmacogenetic studies. Am. J. Psychiatry 170 , 207–217 (2013). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 46. Speed, D. & Balding, D. J. Sumher better estimates the snp heritability of complex traits from summary statistics. Nat. Genet. 51 , 277–284 (2019). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 47. Liu, Y. & Xie, J. Cauchy combination test: a powerful test with analytic p-value calculation under arbitrary dependency structures. J. Am. Stat. Assoc. 115 , 393–402 (2019). [ DOI ] [ PMC free article ] [ PubMed ] 48. Bulik-Sullivan, B. et al. An atlas of genetic correlations across human diseases and traits. Nat. Genet. 47 , 1236–1241 (2015). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 49. Li, X. et al. Dynamic incorporation of multiple in silico functional annotations empowers rare variant association analysis of large whole-genome sequencing studies at scale. Nat. Genet. 52 , 969–983 (2020). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 50. Park, J. et al. Integrated platform for multi-scale molecular imaging and phenotyping of the human brain. Science 384 , eadh9979 (2024). [ DOI ] [ PMC free article ] [ PubMed ] 51. Yao, Z. et al. A high-resolution transcriptomic and spatial atlas of cell types in the whole mouse brain. Nature 624 , 317–332 (2023). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 52. Hyun, J. W. et al. Sgpp: spatial Gaussian predictive process models for neuroimaging data. NeuroImage 89 , 70–80 (2014). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 53. Lin, J.-A. et al. Functional-mixed effects models for candidate genetic mapping in imaging genetic studies. Genet. Epidemiol. 38 , 680–691 (2014). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 54. Zhu, H. et al. Fadtts: functional analysis of diffusion tensor tract statistics. NeuroImage 56 , 1412–1425 (2011). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 55. Yuan, Y. et al. Fmem: Functional mixed effects modeling for the analysis of longitudinal white matter tract data. NeuroImage 84 , 753–764 (2014). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 56. Zhu, H., Fan, J. & Kong, L. Spatially varying coefficient model for neuroimaging data with jump discontinuities. J. Am. Stat. Assoc. 109 , 1084–1098 (2014). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 57. Risk, B. B. & Zhu, H. Ace of space: estimating genetic components of high-dimensional imaging data. Biostatistics 22 , 131–147 (2021). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 58. Shan, Y., Huang, C., Li, Y. & Zhu, H. Merging or ensembling: integrative analysis in multiple neuroimaging studies. Biometrics 80 , ujae003 (2024). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 59. Wang, J. & Li, H. Estimation of genetic correlation with summary association statistics. Biometrika 109 , 421–438 (2022). [ Google Scholar ] 60. Dicker, L. H. Variance estimation in high-dimensional linear models. Biometrika 101 , 269–284 (2014). [ Google Scholar ] 61. Evans, L. M. et al. Comparison of methods that use whole genome data to estimate the heritability and genetic architecture of complex traits. Nat. Genet. 50 , 737–745 (2018). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 62. Hou, K. et al. Accurate estimation of snp-heritability from biobank-scale data irrespective of genetic architecture. Nat. Genet. 51 , 1244–1251 (2019). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 63. Berisa, T. & Pickrell, J. K. Approximately independent linkage disequilibrium blocks in human populations. Bioinformatics 32 , 283 (2016). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 64. Bretherton, C. S., Widmann, M., Dymnikov, V. P., Wallace, J. M. & Bladé, I. The effective number of spatial degrees of freedom of a time-varying field. J. Clim. 12 , 1990–2009 (1999). [ Google Scholar ] 65. Han, X., Xu, C. & Prince, J. L. A topology preserving level set method for geometric deformable models. IEEE Trans. Pattern Anal. Mach. Intell. 25 , 755–768 (2003). [ Google Scholar ] 66. Lorense, W. A High Resolution 3d Surface Construction Algorithm (Association for Computing Machinery, 1987). 67. Wang, Y. et al. Brain surface conformal parameterization with the Ricci flow. IEEE Trans. Med. Imaging 31 , 251–264 (2011). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 68. Wang, Y. et al. Multivariate tensor-based morphometry on surfaces: application to mapping ventricular abnormalities in hiv/aids. NeuroImage 49 , 2141–2157 (2010). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 69. Shi, J. et al. Surface fluid registration of conformal representation: application to detect disease burden and genetic influence on hippocampus. NeuroImage 78 , 111–134 (2013). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 70. Yang, J., Lee, S. H., Goddard, M. E. & Visscher, P. M. Gcta: a tool for genome-wide complex trait analysis. Am. J. Hum. Genet. 88 , 76–82 (2011). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 71. Wightman, D. P. et al. A genome-wide association study with 1,126,563 individuals identifies new risk loci for alzheimer’s disease. Nat. Genet. 53 , 1276–1282 (2021). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 72. Trubetskoy, V. et al. Mapping genomic loci implicates genes and synaptic biology in schizophrenia. Nature 604 , 502–508 (2022). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 73. Grove, J. et al. Identification of common genetic risk variants for autism spectrum disorder. Nat. Genet. 51 , 431–444 (2019). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 74. Mullins, N. et al. Genome-wide association study of more than 40,000 bipolar disorder cases provides new insights into the underlying biology. Nat. Genet. 53 , 817–829 (2021). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 75. Wray, N. R. et al. Genome-wide association analyses identify 44 risk variants and refine the genetic architecture of major depression. Nat. Genet. 50 , 668–681 (2018). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 76. Demontis, D. et al. Genome-wide analyses of adhd identify 27 risk loci, refine the genetic architecture and implicate several cognitive domains. Nat. Genet. 55 , 198–208 (2023). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 77. Nagel, M. et al. Meta-analysis of genome-wide association studies for neuroticism in 449,484 individuals identifies novel genetic loci and pathways. Nat. Genet. 50 , 920–927 (2018). [ DOI ] [ PubMed ] [ Google Scholar ] 78. Savage, J. E. et al. Genome-wide association meta-analysis in 269,867 individuals identifies new genetic and functional links to intelligence. Nat. Genet. 50 , 912–919 (2018). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 79. Jiang, Z. Voxel-level summary statistics of hippocampus shape, white matter microstructure, and cortical surface curvature in UK Biobank ( n = 33,324) 10.5281/zenodo.11404333 (2024). 80. Jiang, Z. Highly efficient imaging genetics [repository] 10.5281/zenodo.18049387 (2025). Associated Data This section collects any data citations, data availability statements, or supplementary materials included in this article. Supplementary Materials Supplementary Information (30.7MB, pdf) 41467_2026_69816_MOESM2_ESM.pdf (118.1KB, pdf) Description of Additional Supplementary Files Supplementary Dataset 1-15 (890.8KB, xlsx) Reporting Summary (4MB, pdf) Transparent Peer Review File (19.3KB, pdf) Data Availability Statement The triplets of summary statistics, as well as the LD matrix (including the pre-computed LD scores in the last column in the. ldinfo file) generated in this study, have been deposited in Zenodo under accession code 13787684 79 , 10.5281/zenodo.13787684. The raw structural MRI and diffusion MRI imaging data, as well as genetic data is under restricted access for privacy issues; access can be obtained by application at https://www.ukbiobank.ac.uk . The data used in this study were obtained from UK Biobank application 22783. The colocalization results were extracted from the NHGRI-EBI GWAS catalog, https://www.ebi.ac.uk/gwas/ . The pre-computed LD scores for the Hapmap3 SNPs are at https://console.cloud.google.com/storage/browser/broad-alkesgroup-public-requester-pays/LDSCORE?pageState=(%22StorageObjectListTable%22:(%22f%22:%22%255B%255D%22)) . Links to summary statistics: schizophrenia, https://figshare.com/articles/dataset/scz2022/19426775?file=34517828 ; Alzheimer’s disease, https://ctg.cncr.nl/documents/p1651/PGCALZ2sumstatsExcluding23andMe.txt.gz ; insomnia, https://ctg.cncr.nl/documents/p1651/insomnia_ukb2b_EUR_sumstats_20190311_with_chrX_mac_100.txt.gz ; educational attainment, https://thessgac.com/papers/3 ; cognitive performance, https://thessgac.com/papers/3 ; autism spectrum disorder, https://pgc.unc.edu/for-researchers/download-results/ ; bipolar disorder https://pgc.unc.edu/for-researchers/download-results/ ; major depressive disorder, https://pgc.unc.edu/for-researchers/download-results/ ; attention deficit hyperactivity disorder, https://pgc.unc.edu/for-researchers/download-results/ ; intelligence, https://ctg.cncr.nl/documents/p1651/SavageJansen_IntMeta_sumstats.zip ; neuroticism, https://ctg.cncr.nl/documents/p1651/sumstats_neuroticism_ctg_format.txt.gz ; depression, https://ctg.cncr.nl/documents/p1651/sumstats_depression_ctg_format.txt.gz ; depression affect subcluster, https://ctg.cncr.nl/documents/p1651/sumstats_depressed_affect_ctg_format.txt.gz ; worry subcluster, https://ctg.cncr.nl/documents/p1651/sumstats_worry_ctg_format.txt.gz . All other data supporting the findings of this study are available in the article and the Supplementary Information files. RVGA is a component in the Highly Efficient Imaging Genetics (HEIG) toolbox, which is a comprehensive solution for voxel-level imaging genetic analysis. HEIG 80 (v1.6.1-alpha) is an open-source Python software, https://github.com/Zhiwen-Owen-Jiang/HEIG (citable version https://zenodo.org/records/18049388 ) with tutorial, https://github.com/Zhiwen-Owen-Jiang/HEIG/wiki and example data, 10.5281/zenodo.14214075; The example code for reproducing the main results, https://github.com/Zhiwen-Owen-Jiang/HEIG/tree/pub/misc/code_manuscript/code.md ; LDSC (v1.0.1), https://github.com/bulik/ldsc ; SumHer, https://dougspeed.com/sumher/ ; GCTA (v1.93.2beta), https://yanglab.westlake.edu.cn/software/gcta/#Overview ; PLINK2 (v2.00a3LM), https://www.cog-genomics.org/plink/2.0/ ); Peaks algorithm (v1.0), https://github.com/wnfldchen/peaks . Articles from Nature Communications are provided here courtesy of Nature Publishing Group ACTIONS View on publisher site PDF (2.9 MB) Cite Collections Permalink PERMALINK Copy RESOURCES Similar articles Cited by other articles Links to NCBI Databases Cite Copy Download .nbib .nbib Format: AMA APA MLA NLM Add to Collections Create a new collection Add to an existing collection Name your collection * Choose a collection Unable to load your collection due to an error Please try again Add Cancel Follow NCBI NCBI on X (formerly known as Twitter) NCBI on Facebook NCBI on LinkedIn NCBI on GitHub NCBI RSS feed Connect with NLM NLM on X (formerly known as Twitter) NLM on Facebook NLM on YouTube National Library of Medicine 8600 Rockville Pike Bethesda, MD 20894 Web Policies FOIA HHS Vulnerability Disclosure Help Accessibility Careers NLM NIH HHS USA.gov Back to Top