Mapping disease critical spatially variable gene programs by integrating spatial transcriptomics with human genetics

Mapping disease critical spatially variable gene programs by integrating spatial transcriptomics with human genetics | bioRxiv /* */ /* */ <!-- <!-- /*! * yepnope1.5.4 * (c) WTFPL, GPLv2 */ (function(a,b,c){function d(a){return"[object Function]"==o.call(a)}function e(a){return"string"==typeof a}function f(){}function g(a){return!a||"loaded"==a||"complete"==a||"uninitialized"==a}function h(){var a=p.shift();q=1,a?a.t?m(function(){("c"==a.t?B.injectCss:B.injectJs)(a.s,0,a.a,a.x,a.e,1)},0):(a(),h()):q=0}function i(a,c,d,e,f,i,j){function k(b){if(!o&&g(l.readyState)&&(u.r=o=1,!q&&h(),l.onload=l.onreadystatechange=null,b)){"img"!=a&&m(function(){t.removeChild(l)},50);for(var d in y[c])y[c].hasOwnProperty(d)&&y[c][d].onload()}}var j=j||B.errorTimeout,l=b.createElement(a),o=0,r=0,u={t:d,s:c,e:f,a:i,x:j};1===y[c]&&(r=1,y[c]=[]),"object"==a?l.data=c:(l.src=c,l.type=a),l.width=l.height="0",l.onerror=l.onload=l.onreadystatechange=function(){k.call(this,r)},p.splice(e,0,u),"img"!=a&&(r||2===y[c]?(t.insertBefore(l,s?null:n),m(k,j)):y[c].push(l))}function j(a,b,c,d,f){return q=0,b=b||"j",e(a)?i("c"==b?v:u,a,b,this.i++,c,d,f):(p.splice(this.i++,0,a),1==p.length&&h()),this}function k(){var a=B;return a.loader={load:j,i:0},a}var l=b.documentElement,m=a.setTimeout,n=b.getElementsByTagName("script")[0],o={}.toString,p=[],q=0,r="MozAppearance"in l.style,s=r&&!!b.createRange().compareNode,t=s?l:n.parentNode,l=a.opera&&"[object Opera]"==o.call(a.opera),l=!!b.attachEvent&&!l,u=r?"object":l?"script":"img",v=l?"script":u,w=Array.isArray||function(a){return"[object Array]"==o.call(a)},x=[],y={},z={timeout:function(a,b){return b.length&&(a.timeout=b[0]),a}},A,B;B=function(a){function b(a){var a=a.split("!"),b=x.length,c=a.pop(),d=a.length,c={url:c,origUrl:c,prefixes:a},e,f,g;for(f=0;f<d;f++)g=a[f].split("="),(e=z[g.shift()])&&(c=e(c,g));for(f=0;f<b;f++)c=x[f](c);return c}function g(a,e,f,g,h){var i=b(a),j=i.autoCallback;i.url.split(".").pop().split("?").shift(),i.bypass||(e&&(e=d(e)?e:e[a]||e[g]||e[a.split("/").pop().split("?")[0]]),i.instead?i.instead(a,e,f,g,h):(y[i.url]?i.noexec=!0:y[i.url]=1,f.load(i.url,i.forceCSS||!i.forceJS&&"css"==i.url.split(".").pop().split("?").shift()?"c":c,i.noexec,i.attrs,i.timeout),(d(e)||d(j))&&f.load(function(){k(),e&&e(i.origUrl,h,g),j&&j(i.origUrl,h,g),y[i.url]=2})))}function h(a,b){function c(a,c){if(a){if(e(a))c||(j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}),g(a,j,b,0,h);else if(Object(a)===a)for(n in m=function(){var b=0,c;for(c in a)a.hasOwnProperty(c)&&b++;return b}(),a)a.hasOwnProperty(n)&&(!c&&!--m&&(d(j)?j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}:j[n]=function(a){return function(){var b=[].slice.call(arguments);a&&a.apply(this,b),l()}}(k[n])),g(a[n],j,b,n,h))}else!c&&l()}var h=!!a.test,i=a.load||a.both,j=a.callback||f,k=j,l=a.complete||f,m,n;c(h?a.yep:a.nope,!!i),i&&c(i)}var i,j,l=this.yepnope.loader;if(e(a))g(a,0,l,0);else if(w(a))for(i=0;i (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0];var j=d.createElement(s);var dl=l!='dataLayer'?'&l='+l:'';j.src='//www.googletagmanager.com/gtm.js?id='+i+dl;j.type='text/javascript';j.async=true;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-M677548'); Skip to main content Home About Submit ALERTS / RSS Search for this keyword Advanced Search New Results Mapping disease critical spatially variable gene programs by integrating spatial transcriptomics with human genetics Hanbyul Lee , Haochen Sun , Xuewei Cao , Berke Karaahmet , Zhijian Li , Hans Ulrich-Klein , Mariko Taga , Gao Wang , Philip L. De Jager , David A. Bennett , Luca Pinello , Xin Jin , Rahul Mazumder , Kushal K. Dey doi: https://doi.org/10.1101/2025.09.24.678397 Hanbyul Lee 1 Computational and Systems Biology, Sloan Kettering Institute, Memorial Sloan Kettering Cancer Center , New York, NY, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site For correspondence: deyk{at}mskcc.org leeh16{at}mskcc.org Haochen Sun 1 Computational and Systems Biology, Sloan Kettering Institute, Memorial Sloan Kettering Cancer Center , New York, NY, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site Xuewei Cao 1 Computational and Systems Biology, Sloan Kettering Institute, Memorial Sloan Kettering Cancer Center , New York, NY, USA 2 Center for Statistical Genetics, The Gertrude H. Sergievsky Center, Columbia University , New York, NY, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site Berke Karaahmet 3 Center for Translational & Computational Neuroimmunology, Columbia University , New York, NY, USA 4 Department of Neurology, Columbia University , New York, NY, USA 5 Taub Institute for Research on Alzheimer’s Disease and the Aging Brain, Columbia University , New York, NY, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site Zhijian Li 6 Broad Institute of MIT and Harvard , Cambridge, MA, USA 7 Molecular Pathology Unit, Center for Cancer Research, Massachusetts General Hospital , Boston, MA, USA 8 Department of Pathology , Harvard Medical School, Boston, MA, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site Hans Ulrich-Klein 3 Center for Translational & Computational Neuroimmunology, Columbia University , New York, NY, USA 4 Department of Neurology, Columbia University , New York, NY, USA 5 Taub Institute for Research on Alzheimer’s Disease and the Aging Brain, Columbia University , New York, NY, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site Mariko Taga 3 Center for Translational & Computational Neuroimmunology, Columbia University , New York, NY, USA 4 Department of Neurology, Columbia University , New York, NY, USA 5 Taub Institute for Research on Alzheimer’s Disease and the Aging Brain, Columbia University , New York, NY, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site Gao Wang 2 Center for Statistical Genetics, The Gertrude H. Sergievsky Center, Columbia University , New York, NY, USA 4 Department of Neurology, Columbia University , New York, NY, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site Philip L. De Jager 3 Center for Translational & Computational Neuroimmunology, Columbia University , New York, NY, USA 4 Department of Neurology, Columbia University , New York, NY, USA 5 Taub Institute for Research on Alzheimer’s Disease and the Aging Brain, Columbia University , New York, NY, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site David A. Bennett 9 Rush Alzheimer’s Disease Center, Rush University Medical Center , Chicago, IL, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site Luca Pinello 6 Broad Institute of MIT and Harvard , Cambridge, MA, USA 7 Molecular Pathology Unit, Center for Cancer Research, Massachusetts General Hospital , Boston, MA, USA 8 Department of Pathology , Harvard Medical School, Boston, MA, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site Xin Jin 10 Department of Neuroscience, Dorris Neuroscience Center , Scripps Research, La Jolla, CA, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site Rahul Mazumder 11 Operations Research Center, Massachusetts Institute of Technology , Cambridge, MA, USA 12 Sloan School of Management, Massachusetts Institute of Technology , Cambridge, MA, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site Kushal K. Dey 1 Computational and Systems Biology, Sloan Kettering Institute, Memorial Sloan Kettering Cancer Center , New York, NY, USA 13 Physiology, Biophysics and Systems Biology, Weill Cornell Medicine , New York, NY, USA 14 Gerstner Sloan Kettering Graduate School of Biomedical Sciences , New York, NY, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site For correspondence: deyk{at}mskcc.org leeh16{at}mskcc.org Abstract Full Text Info/History Metrics Supplementary material Preview PDF Abstract Spatial gene expression patterns underlie tissue organization, development, and disease, yet current methods for detecting spatially variable genes (SVGs) lack the flexibility to capture multi-scale structure, ensure robustness across platforms, and integrate with genetic data to assess disease relevance. We present Spacelink, a unified framework that models spatial variability of a gene at both whole-tissue and cell-type resolution using an adaptive mixture of data-driven spatial kernels and summarizes it using an Effective Spatial Variability (ESV) metric. Spacelink achieved up to 3.2x higher detection power over eight existing global SVG and cell-type SVG methods while showing consistently superior FDR control across 34 different simulation settings and also showed superior cross-platform concordance in matched tissue Visium and CosMx datasets. Applied to 3 healthy CosMx human tissues (brain cortex, lymph node, liver), Spacelink revealed that SVGs are highly informative for 113 complex traits and diseases (average N = 340,406). Spacelink showed up to 2.2x higher disease informativeness over competing methods in tissue-relevant complex diseases and traits, conditional on putative non-spatial expression-level confounders. Applied to a mouse organogenesis Stereo-seq atlas (8 developmental stages), Spacelink identified 145 genes with stage-associated ESV within brain independent of mean expression, that are enriched in pathways like Wnt signaling and Rap1 signaling characterizing early and late development, respectively. Integration with in vivo Perturb-seq targeting 35 de novo ASD risk genes revealed that perturbations in excitatory neurons and astrocytes preferentially altered spatially structured downstream gene programs (1.7– 2.2x higher average ESV across stages than other cell types), many of which were enriched for polygenic autism GWAS loci. In neurodegeneration, analysis of 32 Visium dorsolateral prefrontal cortex samples spanning Alzheimer’s disease (AD) pathology stages identified 334 genes with decreasing ESV along amyloid burden (enriched for glycolysis) and 216 genes with decreasing ESV along tau tangle accumulation (enriched for apoptotic pathways). Several AD risk genes ( PKM, CLU, GPI ) showed conserved reductions in spatial variability with AD pathology in both human and 5xFAD mouse, with PKM linking to a colocalized splicing QTL and amyloid burden QTL variant. These results highlight the utility of Spacelink in decoding spatially variable gene programs that connect tissue architecture to disease genetics. Introduction Spatially resolved transcriptomics (SRT) enable measurement of gene expression at cellular or near-cellular resolution while retaining spatial coordinates, thereby providing an unprecedented opportunity to map spatially variable genes (SVGs) that capture localized regulatory programs and tissue organization 1 – 3 . Linking this spatial variability to human disease genetics holds promise for uncovering tissue microenvironments and cellular niches through which genetic risk factors act, providing new insights into disease mechanisms and informing therapeutic interventions 4 – 7 . Despite their potential, current approaches for detecting SVGs face several limitations that restrict their utility for disease genetics. Standard global SVG methods – nnSVG 8 , SpatialDE 9 , SpatialDE2 10 , SPARK 11 - typically assume a single spatial scale or a limited set of fixed kernels, and thus fail to capture the gene-specific multi-scale organization of expression that arises from layered cell type distributions, fine-grained niches, and broad domains within tissues. Recent cell-type SVG methods (STANCE 12 , Celina 13 ) can capture finer changes in spatial variability within a cell type, however, we highlight that they often fail to appropriately correct for spatial colocalization with other cell types, leading to inflated false positives or reduced power. It is also unclear which metric of spatial variability would be most appropriate for comparing genes in a manner that is stable across platforms 14 , experimental conditions 15 , 16 and physiological contexts 17 – 19 —an essential requirement for integrating SVG programs with human genetics to robustly assess disease relevance. Here, we present Spacelink , a unified statistical framework for detecting and prioritizing SVGs at both global tissue and cell-type resolution. Spacelink employs an adaptive multi-kernel model to capture spatial variance across diverse length scales, and its cell-type-specific version introduces a data-driven gating strategy to correct for spatial colocalization, designed to improve the specificity for cell types that are weakly represented in mixed spots relative to more abundant colocalizing cell types. To summarize spatial variability, we define Effective Spatial Variability (ESV) , a metric which integrates variance magnitude of each component kernel and its corresponding spatial scale into a single interpretable score directly suited for genetic analyses. We extensively benchmarked the global and cell-type versions of Spacelink against 8 state-of-the-art global and cell-type SVG methods across 34 diverse simulation settings informed by both real-world spatial datasets and theoretically distinct spatial patterns; Spacelink consistently outperforms existing methods in power, robustness, and reproducibility. Applications to spatial data from multiple healthy tissue atlases, mouse organogenesis, and Alzheimer’s disease neurodegeneration demonstrate Spacelink’s ability to integrate spatial variability with human genetics, revealing disease-relevant programs that are not captured by conventional expression-based analyses. Results Overview of Methods In Spacelink, we propose a framework to identify and prioritize spatially variable genes using spatially resolved transcriptomics data, both at the resolution of the whole tissue and individual cell types, and evaluate their informativeness in human disease, organogenesis, and tissue degeneration. Spacelink employs multiple data-driven spatial kernels to model the data covariance of each gene ( Figure 1a ). Spacelink models normalized and transformed expression values of a gene g at N spatial locations s = ( s 1 , s 2 , …, s N ) as follows. where X denotes the covariates including the intercept and other spot-level metadata if available, represents kernel bandwidth (spatial length scale), denotes the variance component attributed to kernel , and captures residual noise. By combining kernels with different length scales, Spacelink accommodates spatial variability ranging from fine-scale microenvironments to broad tissue domains. The set of spatial kernels can be viewed as a data-driven selection from a potentially infinitely large collection of spatial kernels of varying bandwidths, which we term as a “ Kernel Store ” ( Methods, Supplementary Figure S1 ). To test whether a gene g is spatially variable, Spacelink evaluates the null hypothesis for all l , against the alternative for some l . We estimate variance components using a Non-Negative Least Squares (NNLS) approach 20 , which links the empirical covariance of each gene to the pre-selected spatial kernels. Beyond hypothesis testing, we introduce a novel gene-level metric, Effective Spatial Variability (ESV), which integrates estimated variance components with their corresponding kernel bandwidths to provide an interpretable score for ranking genes by spatial variability and facilitating genetic integration. Download figure Open in new tab Figure 1. Illustration of the Spacelink framework for identifying spatially variable genes (SVGs) at both global and cell-type levels, and its disease-focused applications. (a). Spacelink selects data-adaptive kernels from a “Kernel Store”—a collection of exponential kernels spanning multiple bandwidths—and models gene expression using a linear mixed model, where total spatial variance is represented as a mixture of kernel-specific variance components. At the global level, SVGs are identified by testing whether all variance components are zero, and spatial variability is quantified using the Effective Spatial Variability (ESV) metric. At the cell-type level, Spacelink extends the model for lower-resolution platforms (e.g., 10x Visium) by (i) incorporating estimated cell type proportions as fixed-effect covariates, (ii) projecting kernel random effects to the focal cell type, and (iii) conditioning out the spatial effects of colocalizing cell types via a data-driven gating procedure. Hypothesis testing and ct-ESV calculation follow the same principles as the global model. (b). Illustrative examples of genes exhibiting diverse global and cell-type level spatial expression patterns. From left to right: (i) a global SVG with high spatial bandwidth of variation resulting in high ESV, (ii) a global SVG showing multiple clusters of signal with low spatial bandwidth of signal, thereby resulting in low ESV, (iii) a non-SVG with random, spotty expression, (iv) a global SVG that is also a ct-SVG in cell type 3, but is not a ct-SVG for cell type 1, based on the cell type distribution in Panel a. (c) Spacelink is applied to SRT data from 3 tissues (brain cortex, lymph node, liver) with a large number of diseases and quantitative traits with well-powered Genome-wide Association Studies (GWAS). We evaluate gene programs based on Spacelink hypothesis test, ESV and ct-ESV in terms of disease informativeness in both matched and non-matched diseases and traits, conditional on putative non spatial confounders. (d). Spacelink is applied to spatiotemporal data along mouse organogenesis and Alzheimer’s disease (AD) human and mouse neurodegeneration data to identify genes whose spatial variability is associated with developmental and degenerative processes. This weighting scheme penalizes kernels with minimal spatial structure, prioritizing genes with strong, non-trivial spatial variability ( Figure 1b ). Genes that do not reach statistical significance in the Spacelink hypothesis test are assigned an ESV score of 0. By construction, ESV downweighs trivial, very short-range autocorrelations that may arise from technical artifacts such as transcript spillover or cell segmentation errors, thereby prioritizing genes with robust spatial organization at biologically meaningful length scales. In lower-resolution SRT platforms such as 10x Visium, where each spot aggregates transcripts from multiple cells spanning different types, it is difficult to resolve spatial variability specific to a single cell type. For example, a weakly representative cell type, like microglia in brain cortex, may colocalize with other dominant cell types, like oligodendrocytes within mixed spots, leading to confounded signals. To this end, we propose a cell-type-specific extension of Spacelink that attributes spatial variability to a focal cell type r, while correcting for potential confounding due to spatial colocalization with other cell types ( Figure 1a ). For a specific cell type r, we model the expression data y g as follows. where denotes fixed effect covariates comprising the intercept, estimated cell type proportions in spots as evaluated by RCTD method 21 , as well as other spot-level metadata if available. Π r represents a N × N diagonal matrix of the proportions of cell type r across N spots. represent random effects capturing the projection of the pre-selected kernels for gene g in the global SVG model to the focal cell type r. A key technical challenge in cell-type SVG modeling is in disentangling the spatial random effects of the focal cell type from that of other colocalizing cell types. This colocalization may artificially introduce correlation among the random effects, and not correcting for it may lead to false positives. On the other hand, naively including random effects for all cell types and length scales can lead to an overcrowding of kernels, where the model becomes overparameterized, diluting the power to detect true spatial signals in the focal cell type. To resolve this, the cell-type-specific version of Spacelink employs a data-driven gating procedure to determine a set of cell types, Coloc(r), that are spatially colocalized with the focal cell type r (see Methods for details) and then corrects for the random effects of each of these cell types at the median bandwidth for gene g. We leverage this model to test whether a gene is spatially variable within a cell type r ( for all l , r is fixed) and use Restricted Maximum Likelihood (REML) approach 22 to generate estimates of . Analogous to ESV , we provide cell-type-specific ct-ESV scores, by replacing with (see Methods for additional details). Genes that do not reach statistical significance in the cell-type-specific Spacelink hypothesis test are assigned a ct-ESV score of 0. The ESV and ct-ESV scores are used as primary metrics to link spatial variability of genes to human disease. We also benchmark ESV against other available spatial metrics, such as PropSV 8 , FSV 9 , model LLR (log-likelihood ratio) 8 , 9 and p-values from SVG hypothesis tests, by evaluating their ability to predict disease-associated genes relevant to the underlying tissue. Spacelink offers several key advantages over existing methods for detecting spatially variable genes (SVGs). First, it provides a unified kernel-based framework for identifying SVGs at both the whole-tissue level and in a cell-type-specific manner ( Figure 1a ) — a capability that, outside of STANCE 12 , is largely missing from other tools. Spacelink selects multiple kernels in a data-driven manner for each gene, which (i) reduces reliance on a limited set of pre-defined kernels as used in SPARK 11 and SpatialDE/SpatialDE2 9 , 10 , and (ii) enables accurate modeling of genes with heterogeneous spatial scales due to multi-scale tissue architecture — a limitation in methods employing a single gene-level spatial length scale such as nnSVG 8 . Second, as highlighted below, Spacelink ESV demonstrates superior transferability across datasets and platforms as well as higher disease informativeness over existing spatial metrics, making it better suited for cross-sample comparisons and downstream phenotypic association. Third, the data-driven random effect gating in the cell type version of Spacelink is designed to improve the specificity and sensitivity of the model compared to no cell type colocalization correction (as in Celina 13 ) and full random effect adjustment across all cell types (as in STANCE 12 ) respectively. Fourth, Spacelink extends the ESV framework to the cell-type level via the ct-ESV score, enabling quantification and prioritization of cell type SVGs; Celina and STANCE do not provide a comparable quantitative metric of cell-type-specific spatial variability. Finally, Spacelink is computationally efficient and scalable to tens of thousands of genes and thousands of spatial locations. It uses closed-form updates for NNLS-based variance component estimation and supports parallelization across genes, enabling rapid inference at large scale. We applied Spacelink to a broad spectrum of spatial transcriptomics datasets from human and mouse, including (i) 3 high single-cell resolution CosMx samples (see Data Availability ) and matched lower-resolution 10x Visium samples from brain cortex, lymph node and liver 23 , 24 ( Data Availability ), (ii) mouse organogenesis Stereo-seq data spanning 8 developmental stages (E9.5 to E16.5) 25 , and (iii) 32 ROSMAP 10x Visium human samples spanning different stages of Alzheimer’s disease pathology 26 ( Data Availability ) and 80 10x Visium mouse samples spanning varying degrees of amyloid deposition in 5xFAD mouse model 6 . We assessed the disease informativeness of Spacelink by using ESV and ct-ESV scores from healthy brain cortex, lymph node, and liver to predict gene-level disease prioritizations from PoPS 27 and MAGMA 28 , 29 across 21, 16 and 31 relatively independent ( r g < 0.8) brain, blood and liver-related diseases and traits ( Figure 1c , Supplementary Table S1 ). Subsequently, we expanded this analysis to 113 complex diseases and traits (average N = 340,406). All disease-related evaluations of Spacelink are conditioned on putative non-spatial confounders such as the mean and overall variance (spatial + non-spatial) gene expression and cell-type specificity of expression. For the organogenesis and AD neurodegeneration, we assess variation in gene-level ESV scores against developmental time and pathological progression in AD samples ( Figure 1d ). We have publicly released Spacelink hypothesis test results, together with ESV and ct-ESV scores at tissue and cell type levels across different datasets ( Data Availability ). Spacelink outperforms other SVG methods in simulation experiments We evaluated Spacelink against 6 leading global SVG detection methods – Moran’s I 30 , nnSVG 8 , SpatialDE 9 , SpatialDE2 10 , SPARK 11 and SPARK-X 31 – that demonstrated the highest predictive accuracy across diverse simulation settings in a recent systematic benchmarking analysis 32 . To ensure robust and representative evaluation, we adopted a framework to design 28 spatial simulation datasets, reflective of realistic spatial patterns. These datasets include (i) scDesign3-based simulations parameterized using 10x Visium data from 10 human samples 33 , (ii) covariance-based Gaussian process capturing multi-scale spatial variance 32 , and (iii) custom-designed simulations representing canonical spatial patterns such as gradients, streaks, and hotspots ( Methods, Supplementary Figure S2 ). Each simulation dataset included 50% non-SVGs out of the total gene set, allowing for a balanced evaluation of sensitivity and specificity. The scDesign3 framework allows varying the degree of spatial variability assigned to simulated SVGs through tunable parameters, which we leveraged to evaluate different SVG ranking methods ( Methods ). For hypothesis testing, we evaluated Spacelink using statistical power and false discovery rate (FDR). For gene prioritization, we used area under the precision–recall curve (AUPRC) across different p-value cutoffs and Kendall’s tau rank correlation between ESV scores and the ground-truth scDesign3 rankings. We benchmarked Spacelink ESV and p-value based rankings against other methods based on method-specific log-likelihood ratio (nnSVG and SpatialDE), hypothesis test p-values (SPARK, SPARK-X and SpatialDE2), Moran’s I coefficient, FSV (SpatialDE) and PropSV (nnSVG). Across diverse simulation scenarios, Spacelink consistently achieved high statistical power (mean: 0.792), outperforming other SVG detection methods by 1.07-3.66x on average across the datasets, while also effectively controlling the FDR below the expected 5% across all simulations ( Figure 2a , Supplementary Figure S3a ). Moran’s I showed comparable power to Spacelink in most scenarios but exhibited poor FDR control, yielding inflated false positives across many simulation scenarios. When evaluated using the Kendall’s tau correlation across 10 scDesign3 simulation scenarios, Spacelink ESV score demonstrated the strongest performance, significantly outperforming all spatial metrics across methods ( Figure 2b ). Notably, ESV outperformed comparable metrics such as FSV (SpatialDE) and propSV (nnSVG), achieving 12.7–37.3% higher correlation across datasets. Moreover, while Spacelink’s global SVG test p-value rankings were slightly less accurate than ESV, they still consistently outperformed p-value-based rankings from other methods by 2.1–12.3% in AUPRC (on average) over the datasets ( Figure 2b , Supplementary Figure S3b ). To assess the robustness of Spacelink to spot dropout, data sparsity, and reduced read depth, we performed (i) spot downsampling at varying probabilities as in Chen et al. 34 , (ii) data sparsification, in which gene expression values were set to zero at selected spatial locations with varying proportions, and (iii) binomial thinning of library size, applied with different probabilities of read-depth reduction following the strategy in Dey et al. 35 (see Methods for details). Spacelink showed robust performance in both SVG identification and gene ranking (based on p-value and ESV) at moderate levels of spot downsampling, data sparsification and library size thinning ( Figure 2c , Supplementary Figure S3c ). In contrast, competing methods showed greater sensitivity to these transformations and degraded more rapidly under extreme conditions. Because nnSVG 8 assumes a single spatial length scale per gene, we generated exponential kernel-based simulation data in which each SVG was simulated with exactly one true length scale – a design that aligns directly with nnSVG’s modeling assumptions. We then assessed performance by comparing the estimated scales from Spacelink and nnSVG to the ground truth, computing the absolute difference between the true length scale and (i) the single scale estimated scale by nnSVG or (ii) the Spacelink-estimated scale with the largest variance contribution ( Methods ). Despite the simulation being tailored to nnSVG’s framework, Spacelink consistently recovered length scales closer to the ground truth over nnSVG (pairwise t-test p-value = 5e-22). nnSVG also exhibited instability under certain spatial configurations ( Figure 2d , Supplementary Figure S3d ); this likely reflects limitations of nnSVG’s nearest-neighbor Gaussian process (NNGP) approximation 36 .? By contrast, Spacelink’s multi-kernel framework is able to capture one or more spatial scales that approximate the true length scale more accurately, even in cases where only a single dominant scale is present in the data. Download figure Open in new tab Figure 2. Spacelink outperforms other global and cell type SVG detection and prioritization methods across a broad range of simulation scenarios. (a). Statistical power and false discovery rate (FDR) of Spacelink versus six other global SVG identification methods across 18 simulation settings informed by real-world spatial data and canonical spatial patterns (results for 10 additional settings are shown in Supplementary Figure S3a ). (b). Left: Kendall’s tau correlation of Spacelink ESV with 8 alternative spatial scoring metrics derived from the 6 methods in Panel (a), across 10 scDesign3-based simulations with known ground-truth ranks. Right: Area under the precision–recall curve (AUPRC) of Spacelink p-values compared with p-values from the six global SVG hypothesis tests across 10 scDesign3-based simulations. Asterisks indicate significance from 1-sided pairwise t-tests comparing Spacelink against each method (*** P < 0.001, ** P < 0.01). (c). Statistical power and FDR of Spacelink against global SVG methods under varying levels of spot downsampling and library size thinning ( Methods) . (d). Left: Comparison of length scale estimates from Spacelink and nnSVG against the ground truth in a set of 3,750 simulated SVGs aligned with nnSVG’s framework. Right: An example simulation where Spacelink closely recovers the true length scale whereas nnSVG does not. (e). Distributions of four cell types in 4 different simulation scenarios constructed based on the human dorsolateral prefrontal cortex (DLPFC) dataset to evaluate ct-SVG methods. (f). AUPRC of Spacelink (ct-SVG) based on hypothesis test p-values compared to Celina and STANCE for each of four cell types across four simulation scenarios. (g). Spacelink ct-ESV distributions for the focal cell type versus other cell types across four scenarios. Numerical results are reported in Supplementary Table S2 . In order to compare the performance of Spacelink with Celina 13 and STANCE 12 for detecting cell type level spatially variable genes (ct-SVGs), we adopted a simulation framework previously proposed in Shang et al. 13 . Under this framework, we simulated single-cell and spot-resolution spatial transcriptomics data from (i) human brain dorsolateral prefrontal cortex (DLPFC) 37 and (ii) mouse primary visual cortex 38 . We assume that cells belong to four different cell types, each mapped to one of four spatial domains corresponding to cortical layers. We designed four distinct scenarios (Scenarios I–IV), varying in the degree of cell type mixing and spatial domain overlap ( Figure 2e , Supplementary Figure S4a ; see Methods for details). For each scenario, we first simulated single-cell resolution expression data for 2,304 genes, including 1,536 ct-SVGs and 384 marker genes that were global SVGs but not ct-SVGs, with varying expression levels across cell types or spatial patterns (see Methods for design considerations including the number of genes). Spot-resolution data were then generated by overlaying square grids onto the single-cell resolution data and aggregating expression counts of cells within each grid to obtain spot-level expression. We applied ct-SVG methods to this spot-resolution data using either oracle or RCTD-estimated cell type proportions 21 , with the simulated single-cell data serving as the reference for RCTD. We benchmarked the predictive accuracy of Spacelink against other methods for each cell type and each scenario using AUPRC. In the human DLPFC simulation using RCTD-estimated cell type proportions, Spacelink achieved significantly higher AUPRC values than Celina and STANCE, with maximum improvements of 2.5- and 4.5%, respectively, across different settings ( Figure 2f ). Similar trends were observed for the mouse visual cortex simulation ( Supplementary Figure S4b ) and when using oracle cell type proportions over estimated RCTD proportions ( Supplementary Figure S4c ). Spacelink’s ct-ESV scores effectively prioritized true ct-SVGs over non-ct-SVGs across all scenarios ( Figure 2g , Supplementary Figure S4d ). Note that due to differences in levels of cell type mixing across layers in the simulation design, the cell type level spatial bandwidths and consequently, the ct-ESV values, are not directly comparable across scenarios. Celina and STANCE were not included in this comparison, as they do not provide a cell-type-specific spatial variability metric like ct-ESV. Finally, since Spacelink differs from Celina and STANCE in (i) employing kernel matrices with data-driven spatial bandwidths and (ii) using a data-driven gating strategy to correct for spatial colocalization in weakly represented cell types, we evaluated the relative performance gain contributed by each of these innovations. The primary performance gain was attributable to the data-driven kernel matrices ( Supplementary Figure S4e ), while the gating mechanism provided an additional modest improvement in performance compared to including all random effects from other cell types (as in STANCE) and excluding them entirely (as in Celina) ( Supplementary Figure S4f ). We performed five secondary analyses. First, we assessed if the addition of other stationary and non-stationary kernels provides improvement over the multi-scale exponential kernels used by Spacelink. Across the scDesign3-based simulation datasets in Figure 2a , our proposed version of Spacelink exhibited higher or comparable statistical power than the expanded model incorporating different Matérn (p=1, 2), Gaussian (Matérn; p= ∞) or periodic kernels, while ensuring FDR control ( Supplementary Figure S5a ). Second, we simulated 10 datasets under the null hypothesis to assess the Type I error control of Spacelink. Spacelink, as well as most other SVG methods, consistently maintained Type I error rates below the 0.05 threshold; however, Moran’s I showed weaker Type I error control compared to other methods across the datasets ( Supplementary Figure S5b ). Third, to assess the robustness of Spacelink to imbalanced data, we generated simulations with varying proportions of non-SVGs, ranging from 5% to 95% of the total gene set. Spacelink demonstrated robust performance and outperformed the other methods across all settings ( Supplementary Figure S5c ). Fourth, using the human DLPFC Visium cortex data with annotated cortical layer information 37 , we evaluated the accuracy of domain detection using the top Spacelink ESV genes, compared to the same number of genes prioritized using other spatial metrics. Spacelink achieved comparable or superior accuracy in identifying spatial domains relative to other methods ( Supplementary Figure S5d ). Lastly, we assessed the computational efficiency of Spacelink by comparing the per-gene run time as a function of the number of spatial spots (ranging from a few hundred to several thousand), relative to other global and ct-SVG methods. The Spacelink hypothesis test exhibited comparable speed to other fast global SVG methods. Among ct-SVG methods, Spacelink was faster than STANCE but slower than Celina ( Supplementary Figure S5e ), which is attributable to Spacelink’s use of a more complex model incorporating an expanded set of random effects in select cases to correct for the effect of other colocalizing cell types. The simulation results demonstrate that Spacelink provides accurate, robust, and interpretable detection and prioritization of spatially variable genes across a wide range of spatial patterns and cell type configurations, consistently outperforming existing methods in statistical calibration and spatial signal recovery. Spacelink outperforms other SVG methods in cross-platform consistency across tissues A desirable quality in a good spatial metric is robustness in gene prioritization across different spatial platforms and assays from matched tissue systems. To assess this, we applied Spacelink to single-cell resolution CosMx spatial transcriptomics data and matched spot-resolution 10x Visium data from 3 healthy human tissues 23 , 24 - brain cortex (CosMx: 188,686 cells, 6,278 genes; Visium: 3,639 spots, 33,525 genes), lymph node (CosMx: 1,852,946 cells, 6,175 genes; Visium: 4,035 spots, 9,116 genes), and liver (CosMx: 332,877 cells, 1,000 genes; Visium: 3,200 spots, 36,602 genes) ( Figure 3a , see Data Availability ). These tissues were selected because they are linked to a broad range of diseases and quantitative traits with well-powered genome-wide association studies (GWAS), making them particularly informative for disease-related benchmarking of Spacelink 39 – 41 . To ensure fair cross-platform comparison, we rasterized the CosMx datasets using SEraster 42 and matched the average number of cells per spot to the corresponding Visium data for each tissue (4.69 cells in brain cortex, 44.13 in lymph node, and 11.06 in liver). We observed strong correlations in ESV scores between CosMx and Visium across all three tissues (r=0.77 in brain cortex, 0.60 in liver, and 0.59 in lymph node) ( Figure 3b ). Consistently across platforms, Spacelink selected fewer spatial bandwidths with non-zero weights for genes in brain cortex compared to lymph node and liver ( Figure 3c , Supplementary Figure S6a ). Moreover, genes in the brain cortex exhibited, on average, 3.3- and 5.6-fold larger normalized spatial length scales—scaled by the minimum distance among spots in each tissue—relative to the lymph node and liver, respectively. These patterns likely reflect the simpler layered organization of gene expression in brain cortex, characterized by long-range spatial variability, in contrast to the more heterogeneous architecture of lymphoid and hepatic tissues ( Supplementary Figure S6b ). We highlight two example genes, CNP and EIF5B , with consistent high and low spatial variability in brain cortex across platforms; in contrast, gene SCGB2A2 shows specific spatial variability only in 10x Visium data ( Figure 3d , Supplementary Figure S6c ). Download figure Open in new tab Figure 3. Spacelink demonstrates superior cross-platform consistency in quantifying spatial variability compared to other methods. (a). Spatial annotations of major cell types in CosMx single-cell resolution SRT data from three human tissues—brain cortex, lymph node, and liver—used as the reference for cross-platform comparison. (b). Pearson correlations of Spacelink global ESV scores, calculated on matched gene panels, between CosMx and Visium data across the three tissues. (c). Proportion of genes with different numbers of spatial length scales with non-zero estimated weights in Spacelink whole tissue resolution models applied to CosMx and Visium data in the three tissues. (d). Scatter plot of global ESV values for all genes in the brain cortex CosMx and Visium panels. (e). Cross-platform consistency of cell type–specific SVG (ct-SVG) methods (Spacelink, Celina, STANCE) between CosMx and Visium panels for two matched tissue systems (left: brain cortex, right: lymph node), assessed by Spearman correlations of p-values for each cell type. Scatter plots show correlations across cell types in brain cortex (left) and lymph node (right), comparing Spacelink against Celina and STANCE. (f). Left: Spearman correlations of Spacelink ESV and five alternative metrics between CosMx and Visium data, restricted to tissue-relevant disease/trait genes in cortex (top) and lymph node (bottom). Right: Average Spearman correlations of p-values from Spacelink, STANCE, and Celina across cell types between CosMx and Visium, restricted to tissue-relevant disease/trait genes in the cortex (top) and lymph node (bottom). In both sets of panels, red circles denote values computed using all genes. Other circles correspond to different tissue-relevant complex diseases and traits. Asterisks indicate significance from 1-sided pairwise t-tests (*** P < 0.001). (g). Fractions of the union of the top 500 PoPS-prioritized genes for tissue-matched traits implicated as top N ranked global SVGs based on Spacelink ESV and 5 other spatial gene prioritization metrics applied to CosMx and Visium data from cortex (top) and lymph node (bottom). We consider 4 values of N (100, 200, 500 and 1,000). Numerical results are provided in Supplementary Table S3 . Next, we evaluated the cell-type-specific version of Spacelink against Celina and STANCE by comparing correlations of p-values obtained from Visium and CosMx data. In Visium, we applied each ct-SVG method using cell type proportions estimated via RCTD 21 with CosMx as the reference ( Supplementary Figure S7a ). For the single-cell-resolution CosMx data, methods were applied directly to cells from the target cell type. Spacelink achieved up to 1.38x, 1.23x, and 1.22x higher cross-platform Spearman correlations than Celina and STANCE in the brain cortex, lymph node, and liver, respectively ( Figure 3e , Supplementary Figure S7b ). These findings underscore the robustness of Spacelink in identifying ct-SVGs across platforms. We further evaluated the cross-platform consistency by focusing on genes relevant to relatively independent ( r g < 0.8) complex diseases and traits linked to each tissue ( Supplementary Table S1 ). We also used this consistency measure as a metric to evaluate Spacelink gene prioritization based on ESV and p-value against other spatial metrics. As primary analysis, we benchmarked Spacelink ESV against other metrics quantifying the extent of spatial variability – nnSVG (PropSV), nnSVG (LLR), Moran’s I coefficient, SpatialDE (FSV), SpatialDE (LLR). As secondary analysis, we compared cross-platform consistency of rankings based on Spacelink p-values against SVG hypothesis test p values from 6 other methods. When restricted to the top 500 tissue-relevant PoPS disease genes 27 , Spacelink ESV and other spatial metrics achieved consistently stronger cross-platform correlations— improving by 10%, 30% and 20% in brain, lymph node and liver respectively ( Figure 3f , Supplementary Figure S8a ). Notably, Spacelink ESV scores exhibited significantly (pairwise t-test, P < 0.01) stronger correlations between CosMx and Visium platforms compared to 5, 4 and 4 out of the 5 competing methods when restricted to tissue-relevant disease genes in brain cortex, lymph node and liver respectively. Spacelink’s p-value-based metric showed lower correlation than the ESV score, it still outperformed all six other p-value-based metrics from competing methods (P < 0.05) in all 3 tissues ( Supplementary Figure S8b ). Spacelink (ct-SVG) achieved significantly higher average cross-platform correlations (across cell types) for disease-prioritized genes compared to STANCE and Celina in brain cortex, lymph node, and liver (P < 0.05) ( Figure 3f , Supplementary Figure S8a ). Together, these results highlight that Spacelink consistently delivers stronger cross-platform reproducibility than other methods, particularly for disease-associated genes across diverse tissue contexts. Finally, we systematically compared the CosMx and Visium spatial platforms in terms of disease informativeness of genes matched across panels. Using Spacelink-detected SVGs at varying rank cutoffs (top 100, 200, 500, and 1000 by ESV), we quantified overlap with the union of the top 500 PoPS-prioritized genes for tissue-matched traits. For this comparison, we subsetted the two panels to the same set of genes for a fair comparison. Across all thresholds, CosMx consistently outperformed Visium, yielding on average 24%, 32%, and 26% higher excess overlap in brain cortex, lymph node, and liver ( Figure 3g , Supplementary Figure S8c ). Similar patterns were observed across other SVG detection methods, underscoring the broader advantage of CosMx over Visium for matched genes for downstream disease association studies; however CosMx panels assay a more limited set of genes compared to Visium’s whole-transcriptome coverage. In summary, Spacelink delivers consistently stronger cross-platform reproducibility of spatially variable genes, especially for disease-relevant genes, outperforming existing methods across a diverse set of tissues. Spacelink prioritized genes are uniquely informative for human diseases We assessed the relevance of spatially variable genes (SVGs) and cell-type-specific SVGs (ct-SVGs), identified by Spacelink in CosMx brain cortex, lymph node, and liver datasets, for 113 complex traits and diseases (average N = 340,406) ( Supplementary Table S1 ). First, we constructed Spacelink-based SVG programs defined by ESV (whole-tissue) and ct-ESV (cell-type-specific) scores. These programs were tested for association with two disease-relevant gene sets: (i) top 500 PoPS-scored genes for brain-, blood-, and liver-related traits, and (ii) MAGMA-prioritized genes (z-score > 3) for the same traits ( Supplementary Table S1 ). Because SVG detection methods often prioritize genes with high expression levels or elevated expression variability, we fit logistic regression models associating each ESV program (X) with disease gene set of each trait (Y) while conditioning on putative confounders (Z) such as mean and variance of log-transformed and normalized expression. To enable balanced evaluation across programs, we computed adjusted Odds Ratios (aORs) and probabilistic recall against the disease gene sets ( Methods ). We benchmarked Spacelink ESV and ct-ESV programs against (i) sc-linker probabilistic cell type programs 43 (genes enriched in specific cell types) and (ii) gsMap cell-level gene specificity scores (GSS) 7 . All programs are non-negative by design, and to ensure fair comparison, each was rescaled to the [0,1] interval ( Methods ). Figure 4a shows the aOR and probabilistic recall of programs with positive and significant regression coefficients (P < 0.05) for each trait. In brain cortex and lymph node, Spacelink global ESV programs achieved 2.6x and 2.2x higher aOR on average and 1.5x and 2.3x higher recall on average compared to sc-linker and GSS programs, while ct-ESV programs showed 5.5x and 4.7x higher aORs but lower probabilistic recall compared to other methods. In liver, only a small number of programs reached significance, likely due to the limited number of genes in the dataset (1,000), which constrained detecting strong trends. Results for ct-ESVs remained consistent when additionally adjusting for the mean and variance of corresponding cell-type level normalized expression ( Supplementary Figure S9a ). Replacing PoPS-scored genes with MAGMA-prioritized genes in the regression models resulted in weaker associations across all programs ( Supplementary Figure S9b ). The restricted gene panels in CosMx and other single-cell–resolution platforms such as Xenium limit the application of genome-wide heritability analyses methods. We note that 60% of cell types implicated fewer than 100 genes even at a relaxed cut-off of ct-ESV (>0.2), all of which when mapped to enhancer-gene links relevant to the tissue, as in sc-linker, generated SNP annotations with very small average annotation size (<0.001); this annotation size is deemed unsuitable for genome-wide heritability regression-based approaches such as S-LDSC 44 , 45 . Download figure Open in new tab Figure 4. Spacelink-prioritized genes are uniquely informative for complex diseases and traits. (a). Logistic regression analyses associating each gene program (probabilistic gene score between 0 and 1) with the top 500 PoPS-prioritized disease gene set for each tissue-matched trait, conditioned on putative non-spatial confounders (overall mean and variance of log-normalized expression), across brain cortex, lymph node, and liver. Gene programs from Spacelink, sc-linker, and gsMap GSS were evaluated using log-transformed adjusted odds ratio (aOR) and probabilistic recall against the disease gene set. Details of the aOR and probabilistic recall calculations are provided in Methods . (b). Method-specific joint regression analyses of all nominally significant (p < 0.05) global and cell-type level Spacelink ESV, sc-linker, and GSS programs in each tissue from the marginal analyses in (a). Similar to (a), the regression model is conditioned on overall mean and variance of log-normalized expression. For GSS, we only include programs that have correlation < 0.5, to ensure less collinearity across programs in the model. Relative AUPRC (against baseline AUPRC; see Methods ) was used to evaluate each joint regression model with respect to each tissue-matched trait. Asterisks indicate significance from 1-sided pairwise t-tests (*** P < 0.001, ** P < 0.01). (c). Full joint regression analyses (“Full”) including all gene programs across the 3 methods, Spacelink, sc-linker, and GSS programs, that are nominally significant (p < 0.05) Spacelink, sc-linker, and GSS programs. For GSS, we only include programs that have correlation < 0.5, to ensure less collinearity across programs in the model. We also consider reduced versions of the full model each excluding sets of programs from each method separately. For example, “Full – GSS” excludes all nominally significant GSS programs from the gsMap method. Relative AUPRC (against baseline AUPRC; see Methods ) was used to evaluate each joint regression model with respect to each tissue-matched trait. Asterisks indicate significance from 1-sided pairwise t-tests (*** P < 0.001, ** P < 0.01). (d). Joint regression analyses of all whole-tissue and cell-type level nominally significant (p < 0.05) spatial gene programs corresponding to Spacelink ESV, nnSVG propSV and SpatialDE FSV; significance is assessed using 1-sided pairwise t-tests (*** P < 0.001, ** P < 0.01). (e). Heatmap of disease or trait-specific enrichment (using log-normalized adjusted odds ratio: log(aOR), conditioned on cell type level overall mean and variance of cell-type resolution sctransform expression data) of 48 Spacelink global ESV and cell-type level ct-ESV programs from three tissues across 113 complex diseases and traits. The top 500 PoPS-prioritized genes were used in the regression model for each trait, as suggested in Weeks et al. 130 . Only results for nominally significant (p < 0.05) log aORs from the disease association analyses are shown. (f). Spatial expression plots of six top-ranked ct-SVGs that are also top PoPS-prioritized genes for relevant traits in brain cortex (left), lymph node (middle), and liver (right). Numbers next to each cell type and trait indicate the gene’s ranks within the full panel based on ct-ESV and PoPS scores, respectively. Numerical results are reported in Supplementary Table S4 . Next, we performed a joint regression analysis of all nominally significant (P < 0.05) Spacelink ESV programs, sc-linker programs, and GSS programs in each tissue from the marginal analysis above, conditioning on the same set of putative confounders as in Figure 4a ( Methods ). Because many GSS programs are highly correlated due to expression similarity across neighboring cells, we restricted the joint analysis to those with pairwise correlations < 0.5. For each disease, we implemented a chromosome-based cross-validation scheme by randomly splitting all genes into two equal groups of training and evaluation set 50 times. For each split, we performed a logistic regression model with a ridge penalty on the training set, and evaluated the performance on the held-out evaluation set. The evaluation metric, relative AUPRC, was defined as the ratio of the model’s test-set AUPRC to the baseline AUPRC ( Methods ). Spacelink ESV programs achieved significantly higher relative AUPRC than GSS programs across all tissues (up to 2.15x higher) and outperformed sc-linker programs in every tissue except liver (up to 2.09x higher; pairwise t-test, P < 0.001) ( Figure 4b ). We next compared a full joint regression model—including all significant Spacelink, sc-linker, and (moderately uncorrelated) GSS programs—with three reduced models, each excluding one program set (Full – Spacelink, Full – sc-linker, Full – GSS). Across tissues, removing Spacelink programs led to significantly lower relative AUPRC than the other models (pairwise t-test, P = 1e-03 or lower) ( Figure 4c ). Notably, models combining Spacelink with other programs (Full, Full – sc-linker, Full – GSS as in Figure 4c ) outperformed Spacelink-only model (as in Figure 4b ), especially in brain cortex. This highlights that Spacelink captures complementary disease-relevant information to sc-linker and gsMap and that combining methods can be advantageous. Next, we benchmarked Spacelink against programs derived from other SVG methods. In the marginal regression analysis of individual programs as in Figure 4a , Spacelink global ESV programs achieved higher probabilistic recall (1.1x and 1.3x) for comparable aOR against PropSV (nnSVG) and FSV (SpatialDE) respectively ( Supplementary Figure S9c ). When applied to individual cell types, Spacelink ct-ESV exhibited 2.8x and 8.6x higher aOR compared to analogous cell-type PropSV and FSV scores ( Supplementary Figure S9c ). In the joint regression model comprising all cell type and global programs for each SVG method, Spacelink ESV programs cumulatively achieved significantly higher relative AUPRC in both cortex and lymph node (pairwise t-test, P < 0.001) ( Figure 4d ). When using binary indicators of the top N SVGs (N = 1,000 in cortex and lymph node; N = 500 in liver) defined by p-value, all six alternative methods—and even Spacelink’s own p-value–based programs—showed markedly lower performance than Spacelink ESV programs in cortex and lymph node ( Supplementary Figure S9d, Supplementary Figure S9e ). In liver, results were similar across methods, likely due to the smaller gene panel; analyses with a broader gene set may uncover clearer performance differences. Finally, we examined individual disease-specific enrichment of 48 Spacelink ESV programs (whole tissue and cell type resolution) from the 3 CosMx tissues against 113 complex diseases and traits, using the adjusted odds ratio (aOR) metric (conditional on mean and variance of normalized expression confounders at the whole-tissue level for whole-tissue programs and at the cell-type level for cell type programs) with respect to top 500 PoPS-prioritized gene sets ( Figure 4e ; Methods ). As expected, we observed strong alignment between tissue-resident ct-SVGs and diseases/traits relevant to the corresponding tissue. Ct-ESV programs in brain and lymph node showed 5.93x and 2.02x higher disease informativeness for related traits, compared to other traits ( Figure 4e ). Ct-ESV programs in liver did not show higher disease informativeness for related traits, likely due to the small gene panel. We also observed several notable cross-tissue disease informativeness of ct-SVGs for diseases and traits – such as brain endothelial cells associated with Coronary Artery Disease (CAD) 46 and blood pressure 47 , liver Hep.3 hepatocyte subtype associated with Schizophrenia 48 , lymph node macrophages associated with lung capacity (FEV1FVC) adjusted for smoking 47 . Endothelial cells form the lining of heart and the broader circulatory system, and are known to become dysfunction in CAD 49 , 50 . Previous studies have shown increased prevalence of liver disease among individuals with Schizophrenia 51 , 52 , and spatial programs in Hep,3 hepatocyte subtype may affect shared pathways like ferroptosis between schizophrenia and metabolic dysfunction 53 . Macrophage subtypes have been implicated in lung fibrosis due to their spatial interaction with tissue-resident fibroblasts, promoting fibrotic processes 54 , 55 . We considered a secondary analysis by replacing PoPS with MAGMA-prioritized (Z score > 3) genes in the disease enrichment analysis in Figure 4e ; the resulting enrichment signals were weaker but still showed consistent tissue-specific patterns ( Supplementary Figure S9f ). We highlight several top-ranked ct-SVGs with high PoPS scores for relevant traits ( Figure 4f ). SNAP25 , ranked 10 th in terms of ct-ESV in brain cortex and among top 150 PoPS-scored genes for Schizophrenia 48 , a synaptic gene associated with multiple neurodevelopmental disorders, is primarily expressed in neurons but has also been shown to be regulated in OPCs 56 . RELN , a top PoPS-prioritized gene (rank 107) for Major Depressive Disorder and the 5th strongest ct-ESV gene in inhibitory neurons (Inh.3 subtype), is implicated to affect specific groups of inhibitory neurons and the development of inhibitory synapses 57 . CR2 , a top 100 PoPS gene (rank=94) for Rheumatoid Arthritis (RA) 58 and a top 10 (rank=6) ct-ESV gene in germinal center B cells, is a known B-cell surface complement receptor whose dysregulation contributes to B-cell hyperactivation and autoantibody production 59 . GDF15 , a top-ranked ct-ESV gene in Hepatocyte 1 (Hep.1) cell subtype in liver ana top 500 PoPS gene for Alanine Aminotransferase (ALT) (rank=313), is a tress- and injury-responsive cytokine, which has been found to be associated with plasma ALT concentrations in obese patients with metabolic dysfunction-associated liver disease 60 , 61 . In summary, our results demonstrate that Spacelink-derived SVG and ct-SVG programs provide uniquely informative signals for complex traits and diseases, complementing existing approaches and revealing biologically plausible gene–cell type associations for disease. Spacelink detects temporal dynamics in spatial variability during mouse organogenesis To investigate how spatial variability in gene expression evolves during embryonic development, we applied Spacelink to single-cell resolution Stereo-seq spatiotemporal transcriptomics data 25 spanning eight stages of mouse organogenesis (E9.5 to E16.5) that capture major developmental transitions ( Figure 5a ). To ensure computational feasibility, we used the “Bin50” resolution dataset (as in previous studies) and further rasterized it with SEraster 42 ( Methods ), yielding spatial spots with an average of 13–14 bins across stages. Because embryos increase in size over time, we corrected for potential systematic biases in raw ESV scores by estimating stage-specific effects and subtracting them from each gene’s ESV values prior to regression analysis ( Methods ). In addition to whole-embryo analyses, we also conducted separate analyses of brain and non-brain compartments, with separate rasterization performed independently for each. This design was motivated by two factors: (i) the brain is robustly represented across all eight stages, consistently accounting for 14–26% of embryo size, whereas other tissues typically cover a much smaller fraction (average 5%) and may be absent at certain stages; and (ii) brain-derived Spacelink SVG programs provided the strongest disease-informative signals compared to other tissues ( Figure 4 ), motivating deeper investigation of its developmental spatial architecture. Due to the absence of well-matched dynamic single cell reference atlases matching the exact stages of mouse organogenesis, we did not perform ct-SVG analysis with Spacelink. Download figure Open in new tab Figure 5. Spacelink detects temporal dynamics in spatial variability during mouse organogenesis. (a). Schematic illustrating mouse embryonic Stereo-seq SRT data across 8 stages (E9.5–E16.5) tracking organogenesis; bottom panel shows spots colored by tissue type. (b). Number of SVGs identified by Spacelink at each stage within the whole embryo (left) and the brain region (right), with bars shaded from light to dark green across developmental time. (c). Volcano plots of gene-level regressions associating stage-corrected ESV (adjusted for embryo size; Methods ) with developmental time in the whole embryo (left) and brain (right). Each point represents the regression slope (x-axis) and –log10 p-value (y-axis). Genes with nominally significant (P < 0.05) negative and positive associations are shown in light and dark green, respectively. (d). Pathway enrichment analyses (based on ConsensusPathdb 63 ) of genes with significantly negative (light green) and positive (dark green) associations between spatial variability and developmental progression in the whole embryo (left) and brain region (right). (e). Spatial expression plots of four genes showing (i) negatively associated ESV with developmental time in whole embryo ( Myl9 ), (ii) positively associated ESV with developmental time in whole embryo ( Rab3a ), (iii) negatively associated ESV in brain with developmental time ( Rspo2 ), and (iv) positively associated ESV in brain with developmental time ( Grin2b ). (f). Average stage-corrected ESV within brain regions at each stage among genes with nominally significantly (P < 0.05) negative (light green) and positive (dark green) associations between spatial variability and developmental progression. (g). Examples of genes showing divergent trends in ESV between brain and non-brain across developmental time. Heatmap shows corrected ESV within brain and non-brain at each stage. Rows are shaded by temporal patterns: light green demonstrates negative ESV trend, dark green demonstrates positive ESV trend, grey demonstrates consistently high ESV throughout development. (h). Proportion of differentially expressed (DE) genes in each tissue among genes which are SVGs within non-brain (top) and brain (bottom) at least 6 of the 8 stages. We rank the tissues based on the decreasing values of the proportion. (i). Disease informativeness of genes with negative (circles) or positive (triangles) temporal ESV associations in brain (top) and non-brain (bottom), focusing on four brain-related Mendelian disorders. Informativeness is measured by the fraction of disorder genes with significant ESV trends (X axis) and the odds ratio between gene sets (Y axis). (j). Spatial expression plots of two monogenic autism risk genes, Nrxn1 and Foxp2 , showing different spatial patterns within brain and non-brain. Nrxn1 exhibits positive ESV trend with development in both brain and non-brain, whereas Foxp2 shows positive ESV trend specifically in brain. Numerical results are reported in Supplementary Table S5 . In the whole-embryo analysis, we expectedly observed a progressive increase in the number of SVGs across developmental stages. For example, Spacelink identified 2.05x more SVGs in E16.5 stage compared to E9.5 stage ( Figure 5b ). Similar patterns were observed in non-brain regions, however, notably for brain, the number of SVGs peaked at the E12.5 stage ( Figure 5b , Supplementary Figure S10a) . To quantify these dynamics at the gene level, we fit a regression model relating stage-corrected ESV values against developmental time ( Methods ) and contrasted it with a control model associating overall mean and variance of sc-transform normalized gene expressions 62 across time. The whole-embryo analysis identified 244 and 51 genes exhibiting (i) nominally significantly positive and negative associations between spatial variability and developmental progression (p_ESV 0.2) ( Figure 5c , Supplementary Figure S10b ). Pathway enrichment analysis using ConsensusPathdb 63 showed that negatively ESV associated genes were enriched in muscle contraction pathways, while positively ESV associated genes were linked to later developmental processes such as Rab regulation of trafficking ( Figure 5d ). In brain regions, 104 and 41 genes showed positive and negative associations between spatial variability and development, also without significant expression shifts (p_ESV 0.2) ( Figure 5c , Supplementary Figure S10b ). 53% and 43% of ESV-associated genes in whole embryo and brain were not identified when associating overall (spatial + non-spatial) variance of normalized expression against developmental time, underscoring the distinct informativeness of spatial variability beyond overall expression variability ( Supplementary Figure S10c) . Negatively ESV-associated genes in brain were enriched in Wnt signaling, whereas positively associated genes involved Rap1 signaling ( Figure 5d ). The Wnt signaling pathway is well known to play a critical role in early vertebrate development by controlling the anterior-posterior axis formation and neural plate patterning, subsequently leading to the formation of neural crest 64 – 66 . Wnt signaling is not among the top 5 pathways (based on p-value < 0.05) when associating mean or variance of log-normalized expression in brain instead of ESV against developmental time ( Supplementary Figure S10d ), highlighting the utility of spatial variability metric. Figure 5e highlights representative genes from the aforementioned pathways ( Myl9 : muscle contraction, Rab3a : Rab regulation of trafficking, Rspo2 : Wnt signaling, Grin2b : Rap1 signaling), illustrating temporal increases or decreases in spatial variability in the whole embryo or brain regions. In the brain, we observed that genes with decreasing spatial variability along development exhibited a marked decline in ESV from E13.5, while genes with increasing spatial variability showed an upward trend beginning at E11.5 ( Figure 5f ). This opposing dynamic explains the observed peak in SVG counts at the E12.5 stage, as early and late developmental spatial programs appear to converge at this stage. Next, we examined differences in spatial variability between brain and non-brain regions. Figure 5g illustrates representative genes with divergent spatial variability dynamics over developmental time in the two different regions. In total, 9,481 and 16,566 genes were identified as SVGs in brain and non-brain, respectively, in at least 6 of the 8 developmental stages. 982 genes showed divergent trends in ESV between brain and non-brain across the developmental time ( Methods, Supplementary Table S5 ; these include key developmental transcription factor genes like Dlx6/Dlx5, Dmrt3 , and Hoxd3. Dlx6 and Dlx5 showed increasing ESV across development in brain (r=0.89 and 0.71, respectively) but decreasing ESV in non-brain regions (r=−0.78 for both). Dlx6 and Dlx5 are homeobox genes that are critical to forebrain development and to the specification and differentiation of inhibitory neurons 67 , 68 . Dmrt3 shows a decreasing trend in ESV across development in brain (r=−0.90) but an increasing trend in non-brain (r=0.91); Dmrt3 plays an essential role in spinal cord development by controlling the differentiation of interneuron subtypes involved in locomotor circuits 69 , 70 . The gene Hoxd3 , on the other hand, shows a decreasing trend in ESV across development in non-brain (r=−0.94) but does not show a noticeable trend in brain (r=−0.10). Hoxd3 is a member of the Hox gene family that governs anterior–posterior patterning and maintains cell fate and identity 71 , 72 . We assessed the overlap of Spacelink SVG gene programs (P < 0.05) in brain and non-brain with differentially expressed (DE) genes in each tissue; in both cases, we see high overlaps with DE genes in tissues like spinal cord, dorsal root ganglion, and choroid plexus ( Figure 5h ). We next evaluated the odds ratio of genes with developmental time–associated ESV in brain and non-brain regions against gene sets corresponding to 4 brain-related Mendelian disorders 73 , 74 . Monogenic autism showed a notably stronger excess overlap with genes exhibiting increasing spatial variability in brain compared to those with decreasing variability (OR = 4.05 vs 1.44) ( Figure 5i ). By contrast, overlap with ESV-associated genes in non-brain was comparatively modest (OR = 1.70 and 1.04 for positively and negatively associated genes, respectively). These findings suggest that spatial processes emerging during later stages of brain development (E12.5 and beyond; Figure 5f ) are specifically important to monogenic autism. On the other hand, genes associated with developmental Mendelian disorders, as well as monogenic psychiatric and neurological disorders, showed relatively weak excess overlap with genes displaying temporal trends in spatial variability in both brain and non-brain regions ( Figure 5i ). We highlight Nrxn1 and Foxp2 , both monogenic autism genes, that exhibit increased spatial variability in late developmental stages in the brain while displaying distinct spatial patterns in non-brain regions ( Figure 5j ). Together, these results highlight that Spacelink can uncover dynamic spatial variability programs during organogenesis that are uniquely informative for monogenic autism and other neurodevelopmental disorders. Spacelink can characterize spatiotemporal dynamics of mouse in-vivo perturbation programs Building on our finding that spatial variability during brain development is particularly informative for monogenic autism, we next examined how perturbations of de novo autism risk genes propagate through spatial programs in the developing brain. To this end, we analyzed an in-vivo Perturb-seq data 75 targeting 35 de novo loss-of-function (LoF) risk genes associated with autism spectrum disorder (ASD) and neurodevelopmental delay (ND) 76 , 77 , where a gene was knocked out (CRISPR-ko) at the E12.5 stage, followed by single-cell RNA-seq readout at the P7 stage ( Figure 6a ). Among downstream-altered genes (which we term as “perturbation program” following Geiger-Schuller*, Eraslan* et al. 78 ) identified across perturbations in five brain cell types—excitatory neurons (Exc), inhibitory neurons (Inh), astroglia, microglia, and oligodendrocytes, we observed strong concordance between astroglia and excitatory neurons (OR = 14.6), compared to an average OR of 2.3 among other cell type pairs ( Figure 6b ). This likely reflects either (i) activation of radial glia within the astroglial compartment toward neurogenic fates upon perturbation 79 , 80 , or (ii) perturbation-induced activation of gene programs underlying astrocyte–glutamatergic neuron interactions 81 , 82 . Overall, we observed a moderate to high concordance in downstream-altered genes across cell types among the 35 autism KO genes (average OR = 5.0) ( Supplementary Figure S11) . We then applied Spacelink to map spatial variability for each perturbed gene and its downstream-altered genes using the spatiotemporal Stereo-seq mouse organogenesis data from Figure 5 . Of the 32 perturbed genes present in the mouse organogenesis dataset, 26 were classified as spatially variable in all 8 developmental stages (mean ESV per gene ranging from 0.07 to 0.80; Supplementary Figure S12a ). Downstream-altered genes in astroglia and excitatory neurons exhibited, on average, 1.76x higher ESV values than the perturbed genes, whereas downstream genes in other cell types (inhibitory neurons, microglia, oligodendrocytes) showed similar ESV levels (to their perturbed genes ( Figure 6c , Supplementary Figure S12b ). Control analyses in two external datasets—(i) genome-wide Perturb-seq in K562 cells 83 and (ii) transcription factor overexpression in embryonic stem cells 84 both showed a 40% and 59% lower average ESV in their downstream perturbation programs compared to the in-vivo perturbation programs in astroglia and excitatory neurons ( Figure 6c ). This highlights that perturbation of autism risk genes preferentially amplifies spatially organized downstream gene programs, specifically within astroglia and excitatory neurons. Next, we investigated genes that were differentially affected downstream of a given perturbation exclusively in top 3 cell type pairs exhibiting the strongest odds ratio in Figure 6b : astrocyte-Exc, microglia-Exc, and astrocyte-microglia. Genes altered in both astroglia and excitatory neurons, but not microglia, showed consistently high stage-corrected ESV across developmental stages ( Supplementary Figure S12c ). In contrast, cell type pairs involving microglia exhibited a gradual increase in ESV from E9.5 to E16.5 stages ( Supplementary Figure S12c ), likely reflecting the progressive intrusion of microglia that typically begins around E9.5 and continues till approximately E14.5 stages 85 . Across the 8 developmental stages, downstream-altered genes from in-vivo Perturb-seq assay showed 10-20% stronger correlation between brain-specific and non-brain-specific ESV values compared to all genes; this improvement remained consistent when restricting the downstream-altered genes to astrocytes and excitatory neurons ( Figure 6d ). Substituting ESV with mean and overall variance of log-normalized expression abolished or inverted these trends ( Supplementary Figure S12d ), emphasizing the unique insight provided by spatial variability. Next, we assessed the association of the ESV of each downstream-altered gene against (i) the number of perturbations affecting the given gene, and (ii) the number of cell types in which we see the gene with altered expression; in both cases, ESV increased significantly with these factors ( Figure 6e ). Among different cell types, the strongest increase in ESV with the number of upstream perturbations was observed in inhibitory neurons ( Supplementary Figure S12e ). Download figure Open in new tab Figure 6. Spacelink characterizes spatiotemporal dynamics of downstream programs in mouse in vivo perturbations of autism risk genes. (a). Schematic of in-vivo Perturb-seq assay from Jin et al. 136 targeting 35 de novo loss-of-function risk genes related to autism spectrum disorder (ASD) and neurodevelopmental delay (ND). Genes were knocked out (CRISPR-ko) at E12.5, followed by single-cell RNA-seq at the postnatal day 7 (P7) stage. Cells were broadly classified into five types—excitatory neurons (Exc), inhibitory neurons (Inh), astrocytes (Astro), microglia (Micro), and oligodendrocytes (ODC)—and downstream-altered genes were identified per perturbation within each cell type. (b). Odds ratios of excess overlap between downstream altered “perturbation program” genes across perturbations for each pair of cell types. Magnitude (odds ratio, dot size) and significance (−log 10 (p-value), dot color) of the excess overlap. (c). (Top panel): Effective Spatial Variability (ESV) of perturbed genes derived from three Perturb-seq experiments across mouse organogenesis stages. These perturbation experiments include our focal in-vivo Perturb-seq assay in brain 136 , as well as two external Perturb-seq assays - genome-wide Perturb-seq in K562 cells 83 , and transcription factor overexpression in embryonic stem cells 84 . Each dot corresponds to a perturbed gene in one of these assays. The points are colored by the underlying perturbation experiment. (Bottom panel): Average ESV of the perturbation program genes for each perturbation. For the in-vivo Perturb-seq, we grouped perturbation programs into two groups – genes altered downstream of a perturbation in astroglia and excitatory neurons (colored maroon), and genes altered downstream of a perturbation in microglia, oligodendrocytes and inhibitory neurons (colored pink). For the external datasets, given a single underlying cell type context, we consider the full set of altered genes in the perturbation program for each perturbation. (d). For each of the 8 mouse organogenesis stages, correlation between brains-specific and non-brain specific ESV across genes from three broad classes: union of perturbation programs of genes in excitatory neurons and astroglia, union of perturbation programs of genes in inhibitory neurons, microglia, and oligodendrocyte, and (iii) all genes. (e). ESV of downstream-altered genes (any cell type) across 8 organogenesis stages stratified by (left) number of perturbed genes affecting them and (right) number of affected cell types. Thick lines denote average ESV; shaded areas show confidence bands. Red dashed line indicates the overall trend estimated by fitting a linear regression model. Asterisks denote significance of Pearson correlation (*** P < 0.001). (f). Scatter plots comparing average ESV of downstream-altered genes (perturbation program) per perturbation at each stage for all genes versus genes intersecting with three autism-associated gene sets: (i) top 100 MAGMA-prioritized polygenic genes (top), top 100 PoPS-prioritized polygenic genes (middle), and (iii) 108 monogenic autism risk genes from OMIM 73 by Freund et al. 74 (bottom). Red lines indicate fitted regression line through the scatter plot based on a linear regression model with zero intercept. We report the slope of the regression line. (g). Spatial expression plots of two de novo autism risk genes, Upf3b and Stard9 across mouse organogenesis stages (left panel), showing minimal spatial variability patterns, and two downstream-altered genes for both these genes (right panel), showing distinct spatial variability patterns in brain and non-brain. Numerical results are reported in Supplementary Table S6 . To assess the polygenic and monogenic disease relevance of spatial variability in perturbation programs from de novo autism risk gene perturbations, we calculated the average ESV of program genes at each developmental stage of mouse organogenesis stratified by three classes of autism risk genes - (i) the top 100 MAGMA-prioritized GWAS genes 28 , 29 , (ii) the top 100 PoPS-prioritized genes (integrating GWAS and gene-level functional features) 27 , and (iii) 108 monogenic autism Mendelian disorder genes 74 . Perturbation program genes overlapping the MAGMA and PoPS sets exhibited markedly higher average ESV (1.47x and 1.36x, respectively) compared to all downstream-altered genes, whereas overlap with Mendelian genes showed no such enrichment ( Figure 6f ). Similarly, relative to the top 100 MAGMA and PoPS gene sets overall, perturbation program genes overlapping these disease gene sets showed stronger ESV (2.43x and 1.75x, respectively; Supplementary Figure S13a ). These enrichments diminished for both MAGMA and PoPS disease gene sets when ESV was replaced by the mean of normalized gene expression ( Supplementary Figure S13b) . Replacing the ESV with overall variance (spatial + non-spatial) resulted in diminished enrichments for MAGMA but showed comparable enrichments for PoPS ( Supplementary Figure S13c ). Together, these findings suggest that perturbing high-effect de novo variants tend to preferably impact spatially variable genes linked to common polygenic autism risk loci. Additionally, perturbation program genes overlapping MAGMA were altered by a greater number of perturbations than those implicated by PoPS or Mendelian autism, with microglial programs showing the strongest associations? ( Supplementary Figure S13d ). In Figure 6g , we highlight two example de-novo autism risk genes, Upf3b and Stard9 , which themselves exhibit low spatial variability in both brain and non-brain regions but influence multiple MAGMA-prioritized autism GWAS genes (Z-score > 3) showing distinct spatiotemporal trends. One such example is Mdk , a downstream target of Upf3b and Stard9 in inhibitory neurons and oligodendrocytes, which shows high spatial variability in early development in brain and non-brain. This is consistent with prior studies demonstrating that Midkines ( Mdk) are strongly expressed in early emryonic development where it assists in neural development, but shows decline in expression in later development stages and in adult tissues 86 , 87 . Another example, Cspg5 , a downstream target of Upf3b and Stard9 in microglia and oligodendrocytes, exhibits increased spatial variability in later stages. We conclude that integrating in-vivo Perturb-seq with Spacelink analysis of spatiotemporal transcriptomic maps enables tracing how genetic disruptions potentially propagate spatially across developmental stages, uncovering convergent spatial programs underlying neurodevelopmental disorders. Spacelink characterizes spatially variable genes programs underlying neurodegeneration To extend our analysis from development to aging and neurodegeneration, we applied Spacelink to 32 ROSMAP 10x Visium samples from brain dorsolateral prefrontal cortex (DLPFC) (4 control samples from 2 healthy donors and 28 samples from 15 donors at different stages of Alzheimer’s disease (AD)) (see Data Availability ), each with detailed global pathological annotations (N =1,166– 2,649 spots, 36,588 genes) ( Figure 7a ). Of the pathological variables, we focused on 3 hallmark indicators of AD progression - the total amyloid plaque burden 88 (square-root transformed), neurofibrillary tangles accumulation (tau) 89 (square-root transformed), and the number of neuritic plaques (square-root transformed). We applied Spacelink to each sample at both whole tissue and cortical layer-level. We did not perform Spacelink cell-type-specific analyses, as disease progression is accompanied by substantial heterogeneity within cell types—such as the emergence of disease-associated microglia 90 , 91 —rendering cell type annotations from static brain cortex snRNA-seq reference-based deconvolution methods like RCTD 21 less reliable. Unlike mouse organogenesis data, the spatial cross sections across samples had similar sizes, leading to consistent length scale initialization ( Supplementary Figure S14a) , thereby precluding the need to correct ESV scores for size differences. Download figure Open in new tab Figure 7. Spacelink captures dynamic alterations in spatial variability across Alzheimer’s disease progression. (a). Illustration of 10x Visium human brain dorsolateral prefrontal cortex (DLPFC) data from 4 control samples from 2 healthy donors and 28 samples from 15 donors at different stages of Alzheimer’s disease (AD), with comprehensive pathological metadata available for each sample, including total amyloid burden, neurofibrillary tangle accumulation, and neuritic plaque counts. Amyloid burden, tangle accumulation, and neuritic plaque counts were square-root transformed to facilitate regression analyses. (b). Boxplots of average Effective Spatial Variability (ESV) scores across genes at control and 3 different stages of AD pathology, determined by levels of amyloid burden (left) and tangle accumulation (right). For each measure, CTRL includes 4 control samples from 2 healthy donors; AD-Lower includes 10 samples with the lowest values of the pathology indicator; AD-Mid includes 9 samples with intermediate values of the pathology indicator; and AD-Higher includes 9 samples with the highest values of the pathology indicator. (c). Volcano plots representing regression analyses associating ESV with amyloid burden (left) or tangle accumulation (right) at the gene level, adjusting for cortical layer composition and mean expression levels, and donor-level random effects. Each point represents the slope coefficient from the regression model and its –log p-value for a gene. We highlight genes that are nominally significant (P < 0.05) for association against amyloid (colored blue) and tau burden (colored orange) respectively. (d). Pathway enrichment analyses of 334 and 216 genes with nominally significant (P < 0.05) negative associations between spatial variability (ESV) and amyloid burden (left) or tangle accumulation (right). Pathway enrichment analysis was performed using ConsensusPathdb 63 with the set of genes in the Visium panel as the background set. (e). Spatial expression plots of the PKM gene (left column), along with layer annotations (middle column), and plaque annotations (right column) in one control sample (top panel) and one high–amyloid-burden sample (bottom panel). (f). Pathology-xQTL colocalization analysis using ColocBoost identifies a 95% colocalized confidence set (chr15: 70531387-70531393) consisting of 5 variants, showing colocalized effect between a splicing QTL phenotype in brain DLPFC and a square-root transformed amyloid burden QTL in an independent ROSMAP cohort with N=595 individuals. (g). Excess overlap analysis (see Methods ) of genes with negative associations between ESV and amyloid burden (from panel c), among genes with the largest average ESV differences between WT and 5xFAD mouse in 10x Visium SRT data from the 5xFAD mouse model at four time points with progressive amyloid accumulation. Asterisks denote significant enrichment by Fisher’s exact test (*** BH-adjusted P < 0.001, ** BH-adjusted P < 0.01, * BH-adjusted P < 0.05). (h) Spatial expression plots of the CLU gene (left column), along with layer annotations (middle column), and plaque annotations (right column) in one control sample (top panel) and one high–amyloid-burden sample (bottom panel) from 10x Visium human DLPFC data. (i) Spatial expression plots of the Clu gene (left column), along with layer annotations (right column) from 10x Visium mouse SRT data of 4-month-old (top panel) and 12-month-old (bottom panel) 5xFAD mice. (j). Alzheimer’s disease (AD)-xQTL colocalization analysis with CLU, showing a 95% colocalized confidence set (chr8: 27604964-27610304) consisting of 8 variants between a CLU eQTL specifically in astrocytes, together with bulk eQTL in DLPFC, and AD GWAS. Numerical results are reported in Supplementary Table S7 . Spacelink revealed a consistent decline in both the number of spatially variable genes (SVGs) and their average ESV scores with increasing amyloid burden, neurofibrillary tangle accumulation, and neuritic plaque counts, indicating a progressive loss of spatial transcriptional organization during neurodegeneration ( Figure 7b , Supplementary Figure S14b ). Across genes, average ESV scores were negatively associated with all AD pathology indicators, even after adjusting for donor-level random effects ( Methods ; Supplementary Figure S14b ). To systematically identify genes whose spatial variability is linked to AD pathology, we performed regression analyses relating gene-level ESVs to amyloid and tau burdens and neuritic plaque counts, while adjusting for cortical layer composition and mean expression levels ( Methods ). We identified 334, 216, and 115 genes whose ESVs were significantly negatively associated with amyloid and tau burdens and neuritic plaque counts, respectively (p < 0.05; Figure 7c , Supplementary Figure S14c ). Pathway enrichment analysis of amyloid-associated SVGs revealed glycolysis and gluconeogenesis as the most enriched pathways ( Figure 7d ), consistent with metabolic dysfunction in early AD. In contrast, genes associated with tangle accumulation were enriched for apoptotic pathways, highlighting neuronal loss and programmed cell death in later disease stages ( Figure 7d ). Genes associated with neuritic plaque counts were enriched for metal ion homeostasis pathways that are known to manage amyloid plaque toxicity 92 ( Supplementary Figure S14d ). We illustrate spatial patterns for two representative genes: PKM , associated with amyloid burden, YWHAG , associated with tau accumulation, and SNCB , associated with neuritic plaque counts ( Figure 7e , Supplementary Figure S14e, Supplementary Figure S14f ). PKM , a key glycolysis enzyme 93 , showed colocalization between a splicing QTL in brain DLPFC and an amyloid burden QTL in an independent ROSMAP cohort of 595 individuals 94 ( Figure 7f ). In contrast to our cross-gene disease association regression model in Figure 4 , we did not include overall variance as a covariate in our primary cross-sample pathology regression. Whereas Figure 4 focused on assessing the unique disease informativeness of entire spatial programs, here our goal was to identify individual genes whose spatial variability is associated with AD pathology. Because overall variance reflects both spatial and non-spatial factors, including it as a covariate risk overcorrecting true SVG signals and reducing power. Nonetheless, we performed a secondary analysis that incorporated overall variance as an additional covariate. The resulting model identified 260, 166, and 103 genes as ESV-associated with amyloid and tau burdens and neuritic plaque counts, respectively, which were enriched in glycolysis, apoptotic, and metal ion homeostasis pathways ( Supplementary Figure S14g, Supplementary Figure S14h ). Finally, in comparing with plaque-proximal differential expression reported by Karaahmet et al. 26 , 25.6% (n=10) of their 39 genes showed significant associations between ESV and at least one of the 3 pathological variables, demonstrating a strong excess overlap between the two classes (19.1x; Methods ). However, the finding that a large fraction (74.4%) of plaque-proximal genes did not exhibit pathology-associated spatial variability underscores that these two spatial gene classes capture complementary aspects of AD-related pathology—one reflecting broad spatial variability and the other localized responses to plaques. Next, we examined mean expression changes associated with AD pathology while adjusting for cortical layer composition ( Methods ). We identified 1,051 and 745 genes whose mean expression levels were significantly negatively associated with tau burden and neuritic plaque counts, respectively (p < 0.05). Only 2 genes showed significant association with amyloid burden ( Supplementary Figure S15a ). The significant genes were enriched in multiple broad neurodegeneration-related pathways, which were distinct from those enriched by ESV-associated genes ( Supplementary Figure S15b ). Next, we applied Spacelink to each of the 6 cortical layers (L1-L6) and white matter (WM) separately, and tested association between layer-specific average ESV and AD pathology markers. We observed nominally significant (p < 0.05) reduction in average ESV with tau accumulation in L6 cortical neurons ( Supplementary Figure S15c ). This is consistent with the fact that tau pathology predominantly affects deep cortical layers such as L6, which are more vulnerable in AD 95 , 96 . However, we observed less pronounced reduction in layer-specific spatial variability compared with global spatial variability for all pathological measures, suggesting that layer-specific spatial architecture is relatively more preserved during AD progression ( Figure 7b , Supplementary Figure S14b, Supplementary Figure S15c ). To evaluate the cross-species robustness of spatial variability patterns identified in aging human brains with AD, we analyzed 10x Visium spatial transcriptomics data from 5xFAD transgenic mouse model 6 , an established model of early-onset familial AD characterized by progressive amyloid plaque deposition 97 . The dataset encompasses 80 samples spanning four time points (4, 6, 8, and 12 months), enabling assessment of spatial variability dynamics along the trajectory of amyloid accumulation. We applied Spacelink to 5xFAD samples from all time points, as well as wild-type (WT) samples, and evaluated the spatial variability patterns of the 334 human amyloid-associated genes observed in Figure 7c . We restricted the Spacelink analysis in mice to the matched cortical layers as AD human brain, and examined differences in average ESV between 5xFAD and WT samples at each time point. We observed higher enrichment of human amyloid-associated ESV genes among the genes with the highest average ESV differences (WT – 5xFAD) at 12 months, the final disease stage, compared to earlier stages with lower amyloid deposition ( Figure 7g ). This is consistent with the human data, where these 334 genes also exhibited reduced spatial variability with increasing amyloid burden. One example is CLU , a known AD risk gene 98 , 99 , which showed 8.48-fold and 1.74-fold lower ESV in the late-stage AD human brain and 12-month-old 5xFAD mouse, respectively, compared to control human samples and early stage mouse 5xFAD samples (4 months) ( Figure 7h , Figure 7i ). Notably, CLU is enriched in expression in astrocytes, and we observed a colocalization between CLU eQTL and AD GWAS specifically in astrocyte cell type 94 ( Figure 7j ); astrocytes are known to play an important role in accumulation of glucose in the vicinity of the amyloid plaques. Other AD risk genes showing conserved spatial variability differences based on amyloid burden include PGK1 and GPI ( Supplementary Figure S15d ). PGK1 is known to support neuronal energy metabolism, and GPI is known to act as both a glycolytic enzyme and a neurotrophic factor 100 , 101 . Some AD risk genes showed distinct ESV trends against amyloid deposition in human and mouse, likely highlighting divergent spatial programs between human AD and mouse 5xFAD models. An example gene is MEF2C which showed a decline in spatial variability in human data but an increasing trend in the 5xFAD mouse model ( Supplementary Figure S15e ). We conclude that Spacelink effectively captures disruption in spatial layered organization of gene expression linked to key pathological features of Alzheimer’s disease—including amyloid and tau burden—and reveals evolutionarily conserved spatial programs between human and mouse models of AD. Discussion We present Spacelink, an adaptive multi-kernel framework for detecting and prioritizing spatially variable genes (SVGs) at both whole-tissue and cell-type resolution. Spacelink introduces three key methodological advances over existing approaches: (i) data-driven estimation of multiple spatial length scales to more accurately model the spatial variance component of each gene; (ii) a rule-based gating strategy for cell-type–specific analysis that corrects for colocalization of weakly represented cell types with other more dominant cell types in mixed spots in lower resolution Visium data, improving both sensitivity and specificity; and (iii) the Effective Spatial Variability (ESV) metric, which integrates spatial length scales and their associated variance components into a single interpretable score. Applied across diverse datasets, Spacelink revealed that spatial variability carries distinctive biological and disease-relevant signals that are not captured by traditional expression-based measures. Spacelink is also able to effectively quantify dynamic changes in spatial variability in organogenesis and tissue degeneration, uncovering regulatory pathways that shift in spatial tissue organization. Finally, Spacelink-derived variability measures provide insights into the preferential spatial activity of gene programs downstream of in vivo perturbations based on cell type context and disease relevance. We note that the ESV metric provides a unified summary of spatial heterogeneity by integrating contributions across multiple spatial scales without making assumptions regarding specific spatial patterns for a gene. In practice, spatial variability in gene expression can arise from diverse processes, including enrichment within spatially defined niches or domains, anatomical gradients, layered structures, and clustered hotspots. By collapsing these heterogeneous signals into a single robust metric, ESV enables consistent prioritization of genes across tissues, datasets, and platforms. Importantly, ESV should not be viewed as an endpoint but more as a prioritization tool: genes with high ESV can subsequently be interrogated with secondary analyses designed to uncover the distinct spatial processes underlying their expression profiles. While previous studies have established SVGs as valuable molecular markers for tissue organization and disease progression 102 , 103 , no consensus has emerged on the optimal metric to quantify spatial variability or its general informativeness across disease contexts. Our results demonstrate that ESV consistently outperforms existing approaches in both robustness and disease relevance, underscoring its utility in connecting spatial tissue architecture to human disease genetics. Our work has several downstream applications. First, Spacelink ESV can be used to pre-filter or weight genes while performing unsupervised clustering or deconvolution of spatial data 104 , 105 and pre-select ligand and receptor pairs for spatial spatial cell–cell interaction inference 106 , 107 , thereby focusing these analyses on genes with strong spatially variable signals. Second, Spacelink tissue-level and cell-type level ESV gene scores can be used as features to inform gene-based disease association statistics 27 , 108 and gene-based probabilistic fine-mapping of transcriptome-wide association studies 109 .Third, Spacelink’s decomposition of spatial variance into multiple length scales enables of genes by their spatial bandwidth profiles across the tissue, which can provide insights into pathway-level spatial organization. Fourth, the relative stability of the Spacelink ESV metric across platforms and conditions demonstrates its utility in comparative analyses across individuals and cohorts, including association with pathology profiles and genetic variation, thereby facilitating spatial quantitative trait loci (QTL) mapping. Our work has several limitations, representing important directions of future research. First, Spacelink prioritizes genes based on overall spatial variability at whole-tissue and cell type resolution, but it is not designed to capture local spatial autocorrelation, as implemented in Hotspot 110 , or model spatial cell-cell interactions 106 , 107 , 111 . As such, Spacelink’s gene prioritization should be viewed as complementary rather than exhaustive in characterizing spatially informative gene programs within a tissue. Second, even with the relative cross-platform stability of ESV over other spatial metrics, residual differences in capture efficiency, spot density, and gene detection rates across platforms and assays can influence downstream comparisons, such as meta-analyses. Third, our analyses primarily focus on platforms profiling the full transcriptome or large gene panels (>1,000 genes), which facilitate genome-wide disease association analysis by providing a sufficiently large set of control genes. To this end, we did not consider highly customized smaller gene panels that are often strongly ascertained based on disease relevance 19 , 112 , 113 . Fourth, although Spacelink decomposes variance into multiple spatial length scales, linking these scales to exact biological structures and mechanisms is non-trivial and may require additional histological support. To this end, pathway enrichment analysis of groups of genes with similar spatial scales can inform about the underlying mechanisms. Fifth, the cell-type–specific version of Spacelink relies on accurate spot-level cell type proportion estimates from deconvolution methods such as RCTD 21 , and errors from poorly matched reference data or cell type annotations can propagate into cell type ESV inference. For this reason, we do not recommend the cell-type version of Spacelink for modeling spatial transcriptomic data from dynamic developmental or degenerative contexts when using static reference datasets, as cell type composition and states may change substantially over time. Sixth, our findings distinguish specific spatially variable programs of genes that play a greater role in genetic risk of disease but do not localize disease risk to a small number of genes, motivating more precise gene-level characterizations. Seventh, similar to other spatial transcriptomic approaches, cell-type-specific version of Spacelink applied to single-cell resolution SRT data will be affected by the cell segmentation quality, in particular, for less prevalent cell types. However, the extent of effect is less pronounced compared to local approaches, and the integration of spatial length scale into ESV score provides an additional safeguard against this artifact. Eighth, while our analyses highlight intriguing synergies between the spatial variability patterns monogenic and polygenic forms of autism, a systematic investigation across multiple polygenic diseases and their matched monogenic counterparts remains to be explored and represents an important direction for future research. Ninth, although we use donor-specific random effects in the AD neurodegeneration data, the way the cortical tissue section has been cut can potentially influence the spatial length scales and ESV, thereby acting as a confounder; While we typically observed robust laminar organization across samples and expect this effect to be modest relative to AD pathology, it remains an important technical consideration when interpreting spatial variability. Tenth, modeling normalized data in Spacelink has limitations, as raw spatial transcriptomics data consist of sparse counts characterized by a high prevalence of zeros and a distinct mean–variance relationship 114 . Directly modeling these large-scale count data, however, is computationally demanding and often prone to convergence issues in high-dimensional settings 115 . Despite all these limitations, Spacelink provides a scalable and interpretable framework for quantifying spatial variability across multiple biological scales, and within and across cell types, enabling robust integration with genetic variation underlying both monogenic and polygenic diseases, and offering a strong foundation for future advances in spatial omics technologies. Methods Overview of Spacelink Spacelink global model We denote y g as the vector of normalized and transformed expression values for gene g across N spatial locations s = ( s 1 , s 2 , …, s N ). At the whole-tissue level, we model y g with a data-driven multi-kernel Gaussian process as follows: Here, X denotes covariates including the intercept and other global metadata at each spatial location; β g the corresponding fixed effects; the spatial variance component for kernel ; and τ g the noise variance. Spacelink assumes that spatially variable genes can exhibit heterogenous spatial patterns with multiple bandwidths due to spatial localization of cell types, spatial zones or domain structure. To capture this, we use a set of exponential covariance kernels , each with a different bandwidth . This set can be viewed as a data-driven subset of a potentially infinite collection of kernels with varying bandwidths, which we term as the “ Kernel Store ”. Theoretically, any completely monotonic isotropic kernel can be approximated by a linear combination of exponential kernels in this store ( Supplementary Note ). We show through secondary simulation experiments ( Supplementary Figure S5a ) that the addition of broader varieties of kernels (Matern, periodic) over exponential kernels do not substantially alter the Spacelink hypothesis test results. Analogous to SpatialDE2 and STANCE, we use a mixture-kernel Gaussian process representation to model gene expression in Spacelink, rather than modeling expression as a mixture of Gaussians with different kernels. This is because the latter is computationally intensive and does not offer substantial advance when assessing global spatial variability, our primary focus, as opposed to local spatial variability. We adopt a 2-stage selection mechanism to adaptively select kernels for each gene. We first pre-select 2 L logarithmically spaced bandwidths ranging from the minimum to the maximum spatial distances across the N locations (with L = 5 as the default) (see schematic in Supplementary Figure S1 ). We then fit the global model for a particular gene g with these candidate bandwidths using Non-Negative Least Squares (NNLS) 20 and identify the smallest and largest bandwidths with non-zero weights (denoted and ). Finally, we define L log-spaced bandwidths within , which specify the final set of spatial bandwidths and corresponding set of kernels . Spacelink assesses whether gene g is spatially variable by testing the null hypothesis for all l , against alternative hypothesis given by for some l . To this end, we first fit the null model using Maximum Likelihood Estimation (MLE) and obtain the score test statistic for each kernel ( Supplementary Note ). Under the null model, each score statistic follows a mixture of chi-square distributions. We apply the Satterthwaite method to approximate the p-value for each statistic and then combine these p-values using the Cauchy combination rule to yield a single, aggregated p-value 116 , 117 ( Supplementary Note ). After obtaining the p-values for all genes, we adjust them using the Benjamini–Hochberg method. The adjusted p-value for each gene is then used to determine whether the gene is spatially variable. Beyond hypothesis testing, we introduce a novel metric— Effective Spatial Variability (ESV) —which integrates variance components with corresponding kernel bandwidths to rank genes by their spatial patterning: The weight penalizes kernels with minimal spatial structure, thereby prioritizing genes with strong, non-trivial spatial variability. Existing metrics such as propSV (σ g / (σ g + τ g ), nnSVG 8 ) and FSV 9 , 32 SpatialDE) lack this property: propSV applies no penalization for spatial bandwidth, while the weighting used in FSV does not effectively filter out genes with minimal spatial structure (e.g., w g = 1 when K g = I in FSV). ESV therefore provides a more interpretable and robust measure of spatial variability. We estimate and by using a Non-Negative Least Squares (NNLS) fit of the empirical covariance matrix: solvable with algorithms such as Lawson-Hanson 20 . Here, is estimated by least squares method. Spacelink cell-type-specific model Spatially variable genes at the whole tissue level are often enriched for cell-type markers due to the non-uniform spatial distribution of different cell types across the tissue. Therefore, it is often desirable to characterize spatial variability not only across the whole tissue, but also within individual cell types. In high-resolution spatial transcriptomics platforms such as CosMx, where each spot corresponds to a single cell with a cell-type annotation, we recommend applying the Spacelink global model to the subset of gene expression restricted to the cells of each cell type for assessing cell-type-specific spatial variability. In contrast, in low-resolution spatial transcriptomics platforms such as 10x Visium, each spot aggregates transcripts from dozens of cells and may contain a mixture of several distinct cell types. In such cases, identifying and quantifying cell type SVGs becomes challenging, as a weakly representative cell type, like microglia in brain cortex, may heavily colocalize with another dominant cell type, like oligodendrocytes, in mixed spots, leading to confounded signals. For such cases, we propose a cell-type-specific extension of Spacelink that models spatial variability attributable to a focal cell type while accounting for potential confounding from spatial colocalization with other celltypes. For a focal cell type r, Spacelink (ctSVG) models the full expression vector y g as a linear mixed model. Here, represents fixed effect covariates comprising the intercept, estimated cell-type proportions in spot as evaluated by the RCTD method 21 as well as other spot-level metadata if available. Π r represents an N × N diagonal matrix of the proportions of cell type r across spots. are the pre-selected kernels for gene g in the Spacelink global model (i.e., the kernels with the L log-spaced bandwidths within ), and denote independent random effects capturing the projection of these kernels onto the cell type r. Coloc(r) denotes a set of cell types that colocalize with the focal cell type r, and our cell-type-specific Spacelink model conditions on the random effect terms for these colocalizing cell types at the median spatial bandwidth m for gene g. A key technical challenge in cell-type SVG modeling is in disentangling spatial random effects of the focal cell type from those of other colocalizing types. This colocalization may artificially induce correlated random effects across cell types and failing to correct for that may lead to false positives. On the other hand, naively including random effects for all cell types and length scales can lead to an overcrowding effect, where the model becomes overparameterized, diluting the power to detect true spatial signals in the focal cell type. To balance these issues, Spacelink (ctSVG) employs a data-driven gating procedure that defines the minimal set of colocalized cell types, Coloc(r), and corrects for their random effects at the median of the pre-selected bandwidths (i.e., the median of the L log-spaced bandwidths within ) for gene g. Restricting correction to the median bandwidth, rather than all bandwidths, avoids excessive removal of true spatial signals from the focal cell type. The set Coloc(r) is determined as follows: Here, denotes the proportion of cell type r at the i -th spot. The first rule states that if the focal cell type r is dominant in most spots where it is present, we do not remove random effects of other cell types. Otherwise, the second rule adaptively identifies the most strongly colocalizing cell types— those more dominant than r in the majority of spots where they are represented—and removes their random effects. We use default cutoffs c 1 = 0.25 and c 2 = 0.2. We considered other cutoffs for c 1 and c 2 in the human DLPFC-based simulation study in Figure 2e , and observed equivalent to worse performance over the default options ( Supplementary Figure S4g ). We assess whether a gene is spatially variable within cell type r by testing the null hypothesis for all l . Under the null, y g follows a normal distribution and we estimate the parameters and τ g via Restricted Maximum Likelihood (REML). Using these estimates, we compute the score test statistic for each kernel, obtain p-values via the Satterthwaite approximation, and then combine them with the Cauchy combination rule. The resulting combined p-value determines whether a gene is spatially variable within cell type r. Analogous to ESV , we suggest cell-type-specific ct-ESV scores, by replacing with , as follows: We estimate and τ g by using Restricted Maximum Likelihood (REML). Simulation design for benchmarking global SVG methods To evaluate global SVG detection methods, we designed three types of simulations: (i) scDesign3-based model 33 , (ii) covariance-based Gaussian process model, and (iii) model using custom-designed canonical spatial patterns. Following recent benchmarking work 32 , we first applied scDesign3 to generate biologically realistic data using 10 human 10x Visium samples as references. For each sample, we selected the top 200 SVGs based on Moran’s I 30 and fitted their expression distributions using the fit_marginal and fit_copula functions in scDesign3. From these, we retained 50 genes with high deviance explained by the fitted model. For each gene, we obtained a mean expression parameter vector reflecting its spatial correlation (denoted μ s ) and generated a non-spatial counterpart ( μ ns ) by randomly shuffling μ s . We then simulated data with varying spatial variability using weighted mixtures α · μ s + (1 − α ) · μ ns as input to the simu_new function, with α ∈ {0.05, 0.1, …,1}. Smaller α values correspond to weaker spatial variability. The α values were used to rank simulated genes, serving as the ground truth for evaluating SVG ranking methods. For each of the 50 genes, we generated 20 replicates across different α values, yielding 1,000 artificial SVGs. An equal number of non-SVGs were generated by shuffling the simulated SVG expressions. Across the 10 Visium datasets, this resulted in 2,000 simulated genes per dataset. Next, we performed covariance-based simulations using a Gaussian process to model multi-scale spatial variance. We considered two versions: (i) a single-kernel covariance and (ii) a mixture of five kernels with different bandwidths. In the single-kernel version, the covariance was defined as while in the five-kernel version, it was defined as Gene expression count data (denoted z ) were simulated using each type of covariance matrix as follows: On the other hand, non-SVGs were simulated separately as described below. Spatial coordinates for all simulations were taken from the 10x Visium Human DLPFC dataset 37 (N = 3,639). The parameter values used in the models were derived from nnSVG applied to the same dataset, which fits the model to normalized and log-transformed gene expression and provides estimates of μ, σ 2 , τ 2 and b for each gene. For the simulation of SVGs, we first selected genes with adjusted p-values < 0.05 according to nnSVG. From these genes, we obtained the mean of μ estimates, the 10th, 30th, 50th, 70th, and 90th percentiles of estimates, and the 10th, 30th, 50th, 70th, and 90th percentiles of τ 2 estimates. These values were used as the and τ 2 parameters in the SVG simulation model. For the single-kernel bandwidth parameter b , among SVGs identified by nnSVG, we used the 25th, 50th, and 75 th percentiles of the b estimates obtained by nnSVG, and generated separate datasets for each b value, which we refer to as ‘1 kernel I’, ‘1 kernel II’, and ‘1 kernel III’, respectively. Forty SVGs were generated for each parameter combination, yielding a total of 1,000 SVGs per dataset. For the five-kernel covariance, three sets of five bandwidth parameter values were used: 10th, 20th, 30th, 40th, and 50th percentiles; 10th, 30th, 50th, 70th, and 90th percentiles; and 50th, 60th, 70th, 80th, and 90th percentiles of b estimates obtained by nnSVG among identified SVGs. These corresponded to the weakest, intermediate, and strongest spatial autocorrelations, respectively, and the resulting datasets were named ‘5 kernels I’, ‘5 kernels II’, and ‘5 kernels III’. Forty SVGs were generated per parameter combination, yielding 1,000 SVGs per dataset. For non-SVGs, the parameter mean of nnSVG estimates for all genes, and τ 2 value was set to the values were set to the 10th, 30th, 50th, 70th, and 90th percentiles of nnSVG estimates among all genes. For each parameter combination, 200 non-SVGs were generated, yielding 1,000 non-SVGs per dataset. Thus, each simulation dataset comprised 2,000 gene expressions. Lastly, we performed custom-designed simulations representing canonical spatial patterns, including gradients, streaks, and hotspots, based on Weber et al. 8 and Shang et al. 13 . As in the covariance-based simulations, spatial coordinates and parameter settings were derived from the 10x Visium Human DLPFC dataset 37 . For non-SVG simulations, we applied the covariance-based design, using the mean of normalized and log-transformed expressions as the parameter μ value, and the 25th, 50th, and 75th percentiles of the variance as τ 2 . For each τ 2 value, 360 non-SVGs were generated, yielding a total of 1,080 non-SVGs per dataset. For SVG simulations, we considered two tissue regions—the white matter and cortical layers, and generated gradient, streak, or hotspot patterns with different bandwidths in each region. For each canonical pattern, we obtained nine distinct spatial patterns having different bandwidths in two regions. After simulating non-SVGs as described above, we modulated expression values along the spatial patterns by downregulating them by 0.25x or 0.5x, or upregulating them by 2x or 4x, thereby generating gene expression data with spatial patterns (see Supplementary Figure S2 ). This procedure produced 1,080 SVGs for each pattern, yielding three datasets with a total of 2,160 simulated gene expression profiles. Additionally, to evaluate method performance under increased noise, we generated versions of each dataset in which 25%, 50%, or 75% of spots were randomly shuffled. To assess the robustness of Spacelink to spot dropout, data sparsity, and reduced read depth, we performed (i) spot downsampling at varying probabilities as in Chen et al. 34 , (ii) data sparsification, in which gene expression values were set to zero at selected spatial locations with varying proportions, and (iii) binomial thinning of library size, applied with different probabilities of read-depth reduction following the strategy in Dey et al. 35 . For all three analyses, we used scDesign3-based simulation data with the 10x Visium Human DLPFC dataset as a reference. In the downsampling experiments, 10%, 20%, …, or 90% of spots were randomly removed across the entire dataset, generating datasets with fewer spatial locations. In the data sparsification experiments, for each gene, 10%, 20%, …, or 90% of spots were randomly selected and their expression values set to zero, producing datasets with varying levels of sparsity. In the binomial thinning experiments, each gene expression count was resampled from a binomial distribution with probabilities of 0.1, 0.2, …, or 0.9. The simulation dataset used to compare length scale estimation between nnSVG 8 and Spacelink was generated using the covariance-based simulation design with a single-kernel covariance. Because both nnSVG and Spacelink operate on normalized and transformed data, we did not use the count data z , but applied the models directly to y . The bandwidth parameter b was set to the 25th, 50th, and 75th percentiles of b estimates among SVGs identified by nnSVG. For each setting, 1,250 SVGs were generated, yielding a total of 3,750 SVGs for the analysis. Simulation design for benchmarking cell type SVG methods To evaluate ct-SVG detection methods, we generated simulated datasets following the framework of Shang et al. 13 , using 10x Visium human DLPFC 37 and STARmap mouse primary visual cortex 38 datasets as references. We first obtained single-cell resolution spatial transcriptomics data and aggregated single-cell measurements into spot-level measurements for the analyses. For the single-cell resolution simulations, we generated 20,000 cells on each tissue using a random-point-pattern Poisson process. Each cell was assigned to one of four spatial domains and one of four cell types, under four scenarios of cell type composition: Scenario 1: Each domain contains only one cell type. Scenario 2: Each domain contains 90% of one cell type and 5% of two others. Scenario 3: Each domain contains 50% of one cell type and 25% of two others. Scenario 4: Each domain contains three cell types at equal proportions. Cell type distributions for each scenario are summarized in Supplementary Table S2 . Based on the original datasets, we simulated gene expression for each cell. A negative binomial model was first fitted to real gene expression count data, and mean and dispersion estimates were extracted. We use the average of mean estimates and the median of the dispersion estimates as the parameter values for the simulation model. Using these parameters, we generated random negative binomial values to simulate non-SVGs (384 per scenario). We next simulated cell type marker genes, defined as genes up- or downregulated only in a specific cell type but without spatial variation. We first generated non-SVGs, and multiplied the expression of cells belonging to the target type by 4, 2, 0.5, or 0.25. For each scenario, 96 marker genes were generated per cell type, totaling 384. To simulate ct-SVGs, we adapted the pattern-based design used for global SVGs, based on a hotspot pattern. For each scenario, we randomly selected one cell of a target cell type as the hotspot center, then defined the hotspot radius using the 25th, 50th, or 75th percentile of distances to other cells of the same type. Within each hotspot, cells were divided into upper and lower halves. After generating either a non-SVG or a marker gene with one of four up- or down-regulation parameter values (4, 2, 0.5, or 0.25), we upscaled expression in the upper half by 2.5 or 2, and in the lower half by 2, or downscaled expression in the upper half by 0.5 or 0.4, and in the lower half by 0.5. This procedure produced ct-SVGs or marker ct-SVGs with hotspot spatial patterns. In each of the 3×4 and 3×4×4 scenarios, 16 ct-SVGs and 4 marker ct-SVGs were generated per cell type, yielding 768 of each in total. As a result, by combining the 384 non-SVGs with 384 non-ct-SVG cell type marker genes, we obtained 768 non-ct-SVGs in total, and for each cell type, 768 ct-SVGs and 768 marker ct-SVGs, thereby allowing for a balanced evaluation of sensitivity and specificity. Finally, to simulate spot-level data, we overlaid a square grid on the single-cell resolution simulation. Each grid cell was treated as a spot, with gene counts summed across all cells inside. From 20,000 cells, we obtained 3,056 spots, each containing between 1 and 19 cells (mean 6–7). We recorded the true cell type proportion of each spot and additionally estimated proportions using RCTD, with the single-cell data as reference. Both the true and inferred proportions were used for ct-SVG methods in evaluation. For simulations under the null hypothesis in one of the secondary analyses, we used 10 scDesign3-based simulation datasets. For each dataset, 1,000 genes were randomly selected, and their expression values were randomly shuffled to generate spatially uncorrelated datasets. To assess robustness to unbalanced data, we used scDesign3-based simulation datasets with the 10x Visium Human DLPFC data 37 as a reference. In one scenario, 1,000 simulated SVGs were retained, and among non-SVGs, 5%, 10%, …, or 95% were randomly retained, producing datasets with more SVGs than non-SVGs. Conversely, 1,000 simulated non-SVGs were retained, and among SVGs, 5%, 10%, …, or 95% were randomly retained, yielding datasets with more non-SVGs than SVGs. For the comparison of computational efficiency, global SVG detection methods were evaluated using the covariance-based simulation design with a single-kernel covariance. Instead of the spatial coordinates from the 10x Visium Human DLPFC dataset 37 , rectangular artificial tissues of various sizes were generated with dimensions 25×20, 25×40, 50×100, and 100×100. For ct-SVG detection methods, scenario 3 of the pattern-based simulation design was used with the 10x Visium Human DLPFC dataset as a reference for parameters. Similarly, rectangular artificial tissues with dimensions 25×20, 25×40, 50×50, 50×100, and 65×100 were used instead of the original coordinates. In each setting, 300 gene expressions were simulated. Performance metrics to benchmark SVG detection methods Throughout the paper, we evaluated SVG detection methods using various metrics, including false discovery rate (FDR), statistical power, area under the precision-recall curve (AUPRC), and Kendall’s tau correlation. FDR and power were computed from p-values with a significance cutoff of 0.05. AUPRC was used to assess ranking performance based on the binary ground truth label (SVG vs. non-SVG), implemented with the pr.curve function from the PRROC package. When a ground truth ranking of spatial variability was available, we used Kendall’s tau correlation for assessment of SVG ranking methods, which is defined as: To assess accuracy in length scale estimation, for nnSVG we computed (true vs. estimated length scale), while for Spacelink we used the following absolute distance: For domain detection methods, we used the adjusted Rand index (ARI) and normalized mutual information (NMI), defined as: where RI is the Rand index measuring the fraction of agreement between two partitions, and where I ( U ; V ) is the mutual information between partitions U and V , and H (·) denotes entropy. Processing spatial transcriptomic datasets and metadata For the analysis of spatial transcriptomic datasets, we used sctransform as the default normalization method. To ensure fairness, all other SVG detection methods were also applied to counts or log-counts normalized with sctransform 62 . Because the focus of our study is not normalization, we did not conduct a systematic comparison of different normalization approaches. It is worth noting that the appropriate strategy for normalizing spatial transcriptomics data—indeed, whether normalization should be performed at all—remains controversial. Prior work has reported that sctransform can sometimes overcorrect for both library size and biological effects, leading to reduced performance in tasks such as spatial domain identification 118 , 119 . However, in our benchmarking against scater 120 , we found that certain technical biases were insufficiently removed by this alternative, resulting in false discoveries in SVG detection. For this reason, we adopted sctransform as the default. Our Spacelink implementation allows users to input data normalized by whichever method they deem most appropriate for their own datasets. In the global SVG simulation analysis, we did not apply library size normalization, as technical biases across spots are absent in simulations. By contrast, in the ct-SVG simulations using spot-resolution data, we applied sctransform because rasterizing cell-resolution data into spots results in a varying number of cells per spot, which introduces variation in library size that requires adjustment. For CosMx datasets and the mouse organogenesis Stereo-seq spatiotemporal dataset 25 , we applied rasterization for computational feasibility using the SEraster::rasterizeGeneExpression function with settings fun = “sum” and square = FALSE 42 . In the CosMx data, we used a resolution of 100 µm, and in the “Bin50” resolution mouse dataset, we used a resolution of 4 for the whole embryo and non-brain regions (rasterization was not applied to brain regions given the smaller number of spots). In the mouse dataset, rasterization was applied without distinguishing between tissues within the whole embryo or non-brain regions, so a single rasterized spot could contain cells from multiple tissues. For analyses comparing CosMx and Visium datasets, we adopted a slightly modified rasterization to better harmonize the two platforms: we used SEraster::rasterizeGeneExpression to obtain spots with neighbors’ distances fixed at 100 µm, then defined a 55 µm-diameter circle around each spot and summed expression values within this circle. This produced rasterized data minimizing discrepancies with Visium. Other SVG methods benchmarked against Spacelink We employed six other methods for global spatially variable gene (SVG) identification and quantification to benchmark the performance against Spacelink. Moran’s I 30 is a standard measure of spatial autocorrelation, where the statistic ranges from –1 to 1; values closer to 1 indicate stronger positive spatial autocorrelation. To compute Moran’s I, we used the spatial_neighbors function from the squidpy 121 Python package with parameters coord_type=“generic” , delaunay=True to define spatial neighbors, followed by the spatial_autocorr function with n_perms=100 to obtain both the Moran’s I statistic and the FDR-corrected one-sided p-value ( pval_z_sim_fdr_bh ). SpatialDE 122 is an early statistical framework for SVG detection that uses Gaussian process regression to model spatial patterns. Using the SpatialDE Python package, we ran SpatialDE.run to obtain the fraction of spatial variance ( FSV ), the log-likelihood ratio ( LLR ), and the corrected p-value ( qval ). SpatialDE2 123 extends the model by combining multiple spatial kernels with a Poisson distribution, providing a framework for spatial count data. It introduces an omnibus test to jointly evaluate multiple kernels under a generalized linear model. We obtained the source code from https://github.com/PMBio/SpatialDE and applied the SpatialDE.test function with omnibus=True to calculate the adjusted p-value ( padj ). nnSVG 124 leverages the nearest neighbor graph to evaluate each gene’s expression pattern across the spatial neighborhood structure, estimating spatial autocorrelation to identify SVGs. Using the nnSVG R package, we ran the nnSVG function with parameters order=“AMMD” , n_neighbors=10 to obtain the proportion of spatial variance ( prop_sv ), the log-likelihood ratio ( LR_stat ), and the adjusted p-value ( padj ). SPARK 125 models count data using a generalized linear spatial model (GLSM) based on the Poisson distribution, testing multiple kernels separately and combining the resulting p-values for SVG detection. SPARK-X 126 adopts a nonparametric covariance test that evaluates multiple kernels with high computational efficiency. We applied the SPARK R package’s spark.test and sparkx functions to obtain the adjusted p-values ( adjusted_pvalue and adjustedPval ) of SPARK and SPARK-X models, respectively. To benchmark performance against the Spacelink cell-type-specific model, we applied two methods—Celina 115 and STANCE 127 —for identifying cell-type-specific spatially variable genes (ct-SVGs). Both rely on Gaussian process models but differ in design: Celina tests multiple kernels and combines the resulting p-values using a Cauchy combination method, while STANCE uses a single kernel. Celina also does not account for the spatial effects of other cell types, whereas STANCE incorporates the spatial effects of all other cell types into the model. For implementation, we followed the official tutorials for each method (Celina: https://lulushang.org/Celina_Tutorial/index.html ; STANCE: https://haroldsu.github.io/STANCE/tutorial.html ).In Celina, we used the Testing_interaction_all function from the CELINA R package to compute the combined p-values ( CombinedPvals ) and determined significance based on an empirical null distribution generated from 10 permutations of spot coordinates, with a 0.01 empirical FDR cutoff, as described in the original paper. In STANCE, we first applied the runTest1 function from the STANCE R package to perform the overall score test. Genes with adjusted p-values ( p_value_adj ) < 0.05 were then subjected to runTest2 to obtain p-values, and significance was assessed at the 0.05 level. We additionally applied two methods to derive gene programs for benchmarking the disease informativeness of the Spacelink gene program. The sc-linker cell type program 43 was generated from scRNA-seq data by performing differential expression analysis for each cell type versus all others using the Wilcoxon rank-sum test, with p-values transformed into scores between 0 and 1. The gsMap gene specificity score (GSS) program 128 was computed by forming micro-domains of homogeneous neighboring spots, averaging gene expression ranks across each micro-domain, and normalizing by the global average rank to yield spot-level specificity scores. All scores were normalized to probabilistic scale ranging from 0 to 1 to facilitate fair comparisons. Processing disease gene-level association statistics To analyze the association of complex traits and diseases with spatial variability, we considered three disease-relevant gene sets: (i) MAGMA-prioritized genes, (ii) PoPS-scored genes, and (iii) Mendelian disorder genes. MAGMA identifies genes with strong GWAS associations by aggregating SNP-level effects into gene-level scores 129 , whereas PoPS trains a linear model on the MAGMA scores using gene-level features derived from single-cell expression data, gene pathways and protein–protein interactions 130 . We used the published PoPS scores as generated by Weeks et al. 130 for 78 UKBiobank traits, and additionally computed the PoPS scores for additional 35 non-UKBiobank traits to facilitate the SVG comparative analysis across tissues ( Supplementary Table S1 ). Mendelian disorder genes for 22 Mendelian disorders, including monogenic autism were curated from OMIM 73 by Freund et al. 74 . For MAGMA, we defined disease-relevant sets by selecting genes with Z-scores greater than 3 or by choosing the top N genes with the highest signed Z-scores, using varying values of N. For PoPS, we defined disease-relevant sets by selecting the top N genes with the highest PoPS scores. Evaluating complex disease informativeness of Spacelink and other SVG programs To assess the disease informativeness of different gene programs, including Spacelink-based SVG programs, while controlling for potential confounders such as mean expression, we fit logistic regression models of the form: Here, y g is a binary variable indicating whether gene g belongs to the top 500 PoPS-scored 27 (or MAGMA-prioritized 28 ) genes for each trait; β 0 is the intercept; X g denotes L predictive features of interest representing gene programs (including Spacelink spatial SVG programs and other programs of interest for benchmarking) with coefficients β 1 ; Z g and contains covariates representing confounders with coefficients β 2 . As confounders, we included the mean and variance of normalized, log-transformed expression values—computed across the whole tissue when comparing different gene programs ( Figures 4a–d ), and within a specific cell type when examining disease-specific enrichment of Spacelink ct-ESV programs ( Figure 4e )—as well as the proportion of non-zero expression values and a binary indicator for high variability (defined by scran::modelGeneVar 131 with FDR < 0.05). We performed two types of analyses. (1) Single-program models Here, X g contained a single gene program ( L = 1) as a predictive feature. Programs included Spacelink ESV or ct-ESV scores, sc-linker programs per cell type 39 , gsMap gene specificity scores (GSS) 7 per cell, nnSVG propSV or ct-propSV scores, SpatialDE FSV or ct-FSV scores, or binary indicators for global or cell-type-specific SVGs defined by SVG detection methods. All these gene programs are non-negative by design and to ensure fair comparison, each program was rescaled to the [0,1] interval as follows. For cell-type-specific programs, SVG detection or scoring methods (including Spacelink) were applied at the cell-type level by subsetting expression in the CosMx data to the corresponding cells. Performance was evaluated using probabilistic recall and adjusted odds ratio, defined as follows: (2) Multi-program models Here, X g contains multiple programs ( L >1) as predictive features in the same model. For example, when evaluating Spacelink scoring method, we constructed an X matrix with the global ESV and all ct-ESV scores (average r=0.57 across datasets); when evaluating sc-linker, we included programs across all cell types (average r=0.18 across datasets). For GSS, due to the large number of programs (equal to the number of spots) and strong correlations among them, we retained only programs with pairwise correlations < 0.5. Genes were partitioned into training and test sets by randomly splitting the 22 autosomes into two groups of 11; genes on one group of chromosomes were assigned to training and the other to testing. We opted for this random split against leave-one-chromosome-out approach as in Weeks et al. 130 and Fabiha et al. 41 , as the number of positives per chromosome is very small. Logistic regression models in this setting were fitted on the training set with a ridge penalty to account for collinearity across features, and performance was assessed on the test set. The evaluation metric was relative AUPRC, defined as the ratio between the model’s AUPRC on the test set and the baseline AUPRC (equal to the proportion of the top 500 PoPS-scored (or MAGMA-prioritized) genes in the test set). Temporal regression of Spacelink ESV in mouse organogenesis spatiotemporal data To investigate the temporal dynamics of gene-level spatial variability while accounting for systematic biases, we employed a two-step regression framework with stage-specific (or sample-specific) random effects. In the first step, we estimated the conditional mean of stage random effect by fitting the following regression model: Here, ESV g,t indicates the raw ESV value for gene g at the t-th developmental stage (t = 1,…,8), α is the intercept, u t represents the random effect at stage t, and ϵ g,t is the noise component. We estimated the conditional mean of random effect, , from this model and subtracted it from each gene’s ESV at each stage: In the second step, we fitted a gene-wise linear regression model using corrected ESV: The same correction and regression procedure was applied to the mean and variance of normalized, log-transformed expression values. Genes were defined as showing an association between spatial variability and developmental progression—but without expression changes—if the slope in the corrected ESV regression was significant (P 0.2). To identify genes with divergent patterns between brain and non-brain re gions, we considered cases where (i) one region exhibited an increasing trend in corrected while the other showed any of the following: a decreasing trend in corrected , a monotonically high profile with corrected ESV > 0.4 in at least 4 stages), or a monotonically low profile ( with corrected ESV < 0.2 in at least 4 stages); or (ii) one region exhibited a decreasing trend while the other showed any of: an increasing trend, a monotonically high profile, or a monotonically low profile, as defined above. Differential ESV analysis in mouse organogenesis spatiotemporal data We conducted differential gene expression (DGE) analysis with 9,481 and 16,566 genes that were identified as SVGs in brain and non-brain regions, respectively, in at least 6 of the 8 developmental stages. For each non-brain tissue, at each developmental stage, we fitted a zero-inflated negative binomial (ZINB) model to test whether genes were expressed at different levels inside the tissue compared to the other non-brain regions. Genes were defined as differentially expressed if they showed a fold change greater than 1 and an adjusted p-value below 0.05 in at least three stages for that tissue. Associating Spacelink against Alzheimer’s disease pathology markers To assess changes in average spatial variability across pathological progression while accounting for donor effects, we fitted the following mixed model: Here, Average ESV d,t denotes the mean ESV value across all genes in the t-th sample from donor d, β 0 and β 1 are fixed effects, x d,t is the pathological measurement for the t-th sample of donor d, u d is the random intercept for donor d, and, ϵ d,t is the residual error term. Using the conditional mean of donor-level random effect , we obtained the donor-corrected average ESV as Average . To systematically identify genes whose spatial variability is associated with Alzheimer’s disease (AD) pathology, we performed regression analyses linking gene-level ESVs across samples to pathological measures, while adjusting for cortical layer composition and mean expression. Specifically, we fit the following model for each gene: where ESV d,t refers to the ESV of the gene in the t-th sample of donor d; y d,t denotes the square root–transformed amyloid burden, tangles, or number of plaques for the sample; Z d,t indicates covariates representing confounders including the mean of normalized and log-transformed expression of the gene and the proportions of spots corresponding to each cortical layer within the sample; u d is the random intercept for donor d. Genes whose regression coefficient for y d,t had a p-value < 0.05 were considered to have spatial variability significantly associated with AD pathology. To identify genes whose mean expression is associated with AD pathology, we applied the same model but replaced ESV g with the mean normalized log-expression of gene g as the predictor, while including cortical layer proportions as covariates. Excess of overlap (EOO) analysis To quantify the concordance between two gene sets, we used an excess of overlap (EOO) statistic 132 . The EOO is defined as the ratio of the observed overlap to the expected overlap under the null hypothesis of random overlap: where n A ∩B is the number of genes shared between the two sets, n A and n B are the sizes of the respective gene sets, and n is the total number of genes considered. We assessed the statistical significance of the excess overlap using Fisher’s exact test. Colocalization analysis We leveraged 17 cis-xQTL contexts spanning 3 molecular modalities (gene expression, splicing, protein abundance) from the aging brain cortex of ROSMAP donors (average N = 595) 133 . All participants enrolled without known dementia and agreed to detailed clinical evaluation and brain donation at death. All studies were approved by an Institutional Review Board of Rush University Medical Center. Each participant signed informed and repository consents and all ROSMAP participants signed an Anatomic Gift Act. The eQTL data encompassed bulk RNA-seq from three cortical regions (dorsolateral prefrontal cortex, posterior cingulate cortex, and head of caudate nucleus), bulk monocytes isolated from peripheral blood, and single-nucleus RNA-seq pseudo-bulk data from six major brain cell types (excitatory neurons, inhibitory neurons, microglia, oligodendrocyte precursor cells, oligodendrocytes, and astrocytes). Bulk-sQTL were generated from three cortical regions using LeafCutter2 134 and bulk-pQTL was performed on the dorsolateral prefrontal cortex. Standard quality control procedures and normalization were applied to all molecular phenotypes, with technical factors including batch effects, RNA integrity, and post-mortem interval adjusted alongside biological covariates such as sex, age at death, and principal components, as described in our previous study 133 . To understand molecular mechanisms underlying spatially variable genes associated with AD and AD pathology, we conducted two targeted multi-trait colocalization analyses using ColocBoost 133 : pathology-xQTL and AD-xQTL colocalizations. For the pathology-xQTL colocalization analysis, we generated pathology phenotypes based on the square-root transformed amyloid burden and tau accumulation metadata corresponding to 595 ROSMAP samples, and then performed the same adjustment by covariates as in cis-xQTL data, except for expression principal components, as the latter can be causally mediated by the underlying pathology. We colocalized these phenotypes against 17 cis-xQTLs in the ColocBoost multi-phenotype colocalization model. For the AD-xQTL colocalization analysis focusing on CLU , we utilized summary statistics from a recent AD GWAS meta-analysis 135 (N_case=111,326, N_control=677,663) after performing quality control against a custom LD reference panel derived from European ancestry samples in the Alzheimer’s Disease Sequencing Project (ADSP) 133 . Data Availability The 10x Visium datasets for DLPFC, prostate, lymph node, heart, kidney, lung cancer, cervical cancer, breast cancer, prostate cancer, and skin melanoma used as input for scDesign3 are available at the following links: https://research.libd.org/spatialLIBD/ (sample ID: 151673), https://www.10xgenomics.com/datasets/normal-human-prostate-ffpe-1-standard-1-3-0 , https://www.10xgenomics.com/datasets/human-lymph-node-1-standard-1-0-0 , https://www.10xgenomics.com/datasets/human-heart-1-standard-1-0-0 , https://www.10xgenomics.com/datasets/human-kidney-11-mm-capture-area-ffpe-2-standard , https://www.10xgenomics.com/datasets/human-lung-cancer-11-mm-capture-area-ffpe-2-standard , https://www.10xgenomics.com/datasets/human-cervical-cancer-1-standard , https://www.10xgenomics.com/datasets/human-breast-cancer-ductal-carcinoma-in-situ-invasive-carcinoma-ffpe-1-standard-1-3-0 , https://www.10xgenomics.com/datasets/human-prostate-cancer-adjacent-normal-section-with-if-staining-ffpe-1-standard , https://www.10xgenomics.com/datasets/human-melanoma-if-stained-ffpe-2-standard , respectively. The STARmap data on mouse visual cortex is available at https://xzhoulab.github.io/SRTsim/ . The CosMx datasets for brain cortex, lymph node and liver are available at the following links: https://nanostring.com/products/cosmx-spatial-molecular-imager/ffpe-dataset/human-frontal-cortex-ffpe-dataset/ , https://nanostring.com/products/cosmx-spatial-molecular-imager/ffpe-dataset/cosmx-human-lymph-node-ffpe-dataset/ , https://nanostring.com/products/cosmx-spatial-molecular-imager/ffpe-dataset/human-liver-rna-ffpe-dataset/ The 10x Visium dataset for lymph node used for cross-platform evaluation is available at https://www.10xgenomics.com/datasets/human-lymph-node-1-standard-1-1-0 . The Stereo-seq mouse organogenesis spatiotemporal transcriptomics data is available at https://db.cngb.org/stomics/mosta/ . The mouse in-vivo Perturb-seq data is available through the Gene Expression Omnibus (accession GSE157977). The 32 ROSMAP 10x Visium samples from brain dorsolateral prefrontal cortex (DLPFC) at different stages of Alzheimer’s disease (AD) are available via the AD Knowledge Portal ( https://adknowledgeportal.org ). Pathologic and phenotypic data from ROSMAP are available at https://www.radc.rush.edu . The 10x Visium data from 5xFAD transgenic mouse model is available through the Gene Expression Omnibus (accession GSE233208). The supplementary tables contain the Spacelink hypothesis test results, ESV and ct-ESV scores across different datasets, as well as the data needed to reproduce the main and supplementary figures. Code Availability Spacelink is available as an R package, with the source code freely accessible at https://github.com/hanbyul-lee/spacelink . Supplementary Figures Supplementary Figure S1. Schematic for selecting data-adaptive length scales . We first initialize 2 L (default L = 5) log-spaced bandwidths between the minimum and maximum spatial distances across N locations, which are shared across all genes. For each gene g, we fit the Spacelink global model using Non-Negative Least Squares (NNLS) and identify the smallest and largest bandwidths with non-zero weights ( and ). We then define L log-spaced bandwidths within to obtain the final gene-specific set of spatial bandwidths and corresponding kernels. Supplementary Figure S2. Schematic of scDesign3-based and pattern-based simulations . (a). scDesign3-based simulation . We first selected the top 50 SVGs using Moran’s I and model deviance explained (via scDesign3 functions). For each gene, we derived a mean expression parameter vector capturing spatial correlation ( μ s ) and generated a non-spatial counterpart ( μ ns ) by randomly shuffling μ s . Data were simulated with varying spatial variability using weighted mixtures α · μ s + 1(1 − α ) · μ ns , with α ∈ {0.05, 0.1, …, 1}. Smaller α values indicate weaker spatial variability, while larger values indicate stronger spatial variability. (b). Pattern-based canonical simulation . We used two tissue regions from 10x Visium Human DLPFC data—the white matter and cortical layers— and generated hotspot, streak, or gradient with different bandwidths. Each canonical pattern (hotspot, streak, gradient) produced nine distinct spatial configurations (three bandwidths × two regions; top panel). By up- or downregulating non-SVG expression along these spatial patterns, we generated gene expression data with spatial structure (bottom panel, left). To assess robustness, we additionally generated versions of each dataset in which 25%, 50%, or 75% of spots were randomly shuffled (bottom panel, right). Supplementary Figure S3. Additional results from global SVG simulation analyses . (a). Statistical power and false discovery rate (FDR) of Spacelink versus six other global SVG identification methods across 6 covariance-based and 12 pattern-based simulation settings. (b). Area under the precision–recall curve (AUPRC) of Spacelink p-values compared with p-values from the six global SVG hypothesis tests across 6 covariance-based and 12 pattern-based simulations. Asterisks indicate significance from 1-sided pairwise t-tests comparing Spacelink against each method (*** P < 0.001). (c). Statistical power, FDR, AUPRC, and Kendall’s tau correlation of Spacelink versus global SVG methods under varying levels of spot downsampling, data sparsification, and library size thinning ( Methods ). (d). An alternative visual illustration of Figure 2d . Left: Violin plots of distances between true and estimated length scales from Spacelink and nnSVG in a set of 3,750 simulated SVGs aligned with nnSVG’s framework. Right: An example simulation where Spacelink closely recovers the true length scale whereas nnSVG does not. Numerical results are reported in Supplementary Table S2 . Supplementary Figure S4. Additional results from cell-type-specific SVG simulation analyses . (a). Distributions of four cell types in 4 different simulation scenarios constructed based on the mouse visual cortex dataset to evaluate ct-SVG methods ( Methods ). (b). Statistical power, false discovery rate (FDR), and area under the precision–recall curve (AUPRC) of Spacelink (ct-SVG) compared with Celina and STANCE for each of the four cell types across four simulation scenarios based on the human dorsolateral prefrontal cortex (DLPFC) dataset (left) and the mouse visual cortex dataset (right). Cell type proportions estimated by RCTD were used to apply each ct-SVG method. (c). Statistical power, FDR, and AUPRC of Spacelink (ct-SVG) compared with Celina and STANCE for each of the four cell types across four simulation scenarios based on the human DLPFC dataset (left) and the mouse visual cortex dataset (right), using oracle cell-type proportions derived from the true cell-type distributions. (d). Spacelink ct-ESV distributions for the focal cell type versus other cell types across four scenarios from the mouse visual cortex dataset. (e). AUPRC of two variants of the Spacelink method versus Celina and STANCE for each of the four cell types across four simulation scenarios from the human DLPFC dataset (top) and the mouse visual cortex dataset (bottom). Spacelink variant 1 and variant 2 use the same data-driven multiple kernels as the original Spacelink (ct-SVG) but differ in how random effects from non-focal cell types are modeled: variant 1 excludes these effects completely (analogous to Celina), whereas variant 2 includes all of them (analogous to STANCE). Asterisks indicate significance from 1-sided pairwise t-tests comparing Spacelink (variant or original) with Celina or STANCE (*** P < 0.001, ** P < 0.01, * P < 0.05). (f). AUPRC of the Spacelink variants described in (e) for each of the four cell types across four simulation scenarios from the human DLPFC dataset (top) and the mouse visual cortex dataset (bottom). (g). Average FDR, statistical power, and AUPRC of Spacelink across four cell types and four simulation scenarios in the human DLPFC dataset, using varying cutoffs c 1 and c 2 to define the set of colocalized cell types, Coloc(r) (see Methods ). The x-axis indicates different values of c 1 and circle color indicates different values of c 2 . The red dashed line marks the default setting (c 2 = 0.25, c 1 = 0.2). Numerical results are reported in Supplementary Table S2 . Supplementary Figure S5. Results of secondary analyses benchmarking Spacelink against other SVG methods . (a). Statistical power and false discovery rate (FDR) of the original Spacelink and 4 variants of the Spacelink method incorporating an expanded kernel set, such as Matérn kernels (p = 1 or 2), Gaussian kernels (equivalent to Matérn with p = ∞), or cosine periodic kernels, across 10 scDesign3-based simulations. (b). Type I error rates of Spacelink compared with six other global SVG identification methods across 10 null simulations generated from scDesign-3 based simulations spanning 10 Visium tissues by randomly shuffling the gene expression across spots. (c). Statistical power, FDR, and AUPRC of Spacelink versus global SVG methods in unbalanced datasets ( Methods ). (d). Adjusted Rand Index (top) and Normalized Mutual Information (bottom) of Spacelink ESV and nine other gene scoring metrics used to select top-ranked genes for domain detection of cortical layers in the human DLPFC Visium datasets 23 . Different boxplot colors indicate the number of top-ranked genes used. Asterisks indicate significance from 1-sided pairwise t-tests comparing Spacelink ESV with each other method (*** P < 0.001, ** P < 0.01, * P < 0.05). (e). Runtime of Spacelink versus six other global SVG identification methods (left) and two other cell-type SVG identification methods (right) for processing expression data of a single gene with varying numbers of spots. Numerical results are reported in Supplementary Table S2 . Supplementary Figure S6. Gene-specific spatial bandwidths estimated by Spacelink in three human tissues, cortical layer structure in the brain cortex, and examples of gene expression in CosMx and Visium brain cortex . (a). Heatmap of the weights (i.e., estimated variance components ) corresponding to gene-specific length scales inferred by the Spacelink model for each gene in brain cortex (left), lymph node (middle), and liver (right) tissues from CosMx (top) and Visium (bottom) datasets. Across spatial platforms, we observe greater heterogeneity in spatial length scales in lymph node and liver, compared to cortex. (b). Cortical layer structure of the brain cortex in the Visium dataset. (c). Spatial expression plots of three genes ( CNP, SCGB2A2, EIF5B ) illustrating distinct ESV patterns between CosMx and Visium in brain cortex (also highlighted in the scatter plot in Figure 3d ). CNP and EIF5B exhibit consistently high and low ESV values across both platforms, whereas SCGB2A2 shows high spatial variability in Visium but not in CosMx. Numerical results are provided in Supplementary Table S3 . Supplementary Figure S7. Estimated cell type proportions in Visium human samples and cross-platform evaluation of ct-SVG methods in liver . (a). Cell type proportions in the brain cortex (top), lymph node (middle), and liver (bottom) Visium datasets, estimated by RCTD using the corresponding CosMx datasets as references. (b). Cross-platform consistency of cell-type-specific SVG (ct-SVG) methods (Spacelink, Celina, STANCE) between CosMx and Visium in the liver, evaluated by Spearman correlations of p-values for each cell type. Scatter plots show correlations across liver cell types, comparing Spacelink with Celina and STANCE. Numerical results are provided in Supplementary Table S3 . Supplementary Figure S8. Additional results from cross-platform evaluation between matched tissue CosMx and Visium panel data . (a). Left: Spearman correlations of Spacelink ESV and five alternative metrics of spatial variability between CosMx and Visium data, restricted to tissue-relevant disease/trait genes in liver. Right: Average Spearman correlations of p-values from Spacelink, STANCE, and Celina across cell types between CosMx and Visium, also restricted to tissue-relevant disease/trait genes in liver (we did not benchmark ct-ESV as Celina and STANCE do not offer a comparable metric). Red circles denote values computed using all genes; other circles correspond to individual tissue-relevant complex diseases and traits. Asterisks indicate significance from 1-sided pairwise t-tests (*** P < 0.001, ** P < 0.01, * P < 0.05). (b). Spearman correlations of p-values from Spacelink and six other SVG methods between CosMx and Visium data, restricted to tissue-relevant disease/trait genes in brain cortex (top left), lymph node (top right), and liver (bottom). Red circles denote values computed using all genes; other circles correspond to individual tissue-relevant complex diseases and traits. Asterisks indicate significance from 1-sided pairwise t-tests (*** P < 0.001, ** P < 0.01, * P < 0.05). (c). Fractions of the union of the top 500 PoPS-prioritized genes for tissue-matched traits implicated as top N ranked global SVGs based on Spacelink ESV and 5 other spatial gene prioritization metrics applied to CosMx and Visium data from liver. We consider 4 values of N (100, 200, 500 and 1,000). Numerical results are provided in Supplementary Table S3 . Supplementary Figure S9. Additional results from disease-informativeness analyses of gene programs in 3 CosMx human tissues – brain cortex, lymph node and liver . (a). Logistic regression analyses as in Figure 4a , but for Spacelink ct-ESV programs, conditioned on additional cell-type-level confounders (mean and variance of log-normalized expression within each cell type). (b). Logistic regression analyses as in Figure 4a , but using MAGMA-prioritized disease gene set (z-score > 3) for each tissue-matched trait, instead of the top 500 PoPS-prioritized disease gene set. (c). Logistic regression analyses as in Figure 4a , but using gene programs obtained from Spacelink ESV, nnSVG propSV, and SpatialDE FSV. (d). Logistic regression analyses as in Figure 4a , but using gene programs obtained from Spacelink ESV and p-values from Spacelink and six other SVG identification methods. (e). Method-specific joint regression analyses of all nominally significant (p < 0.05) global and cell-type-level gene programs in each tissue, as in Figure 4b , but using gene programs obtained from Spacelink ESV and p-values from Spacelink and six other SVG identification methods. Relative AUPRC (against baseline AUPRC; see Methods ) was used to evaluate each joint regression model with respect to each tissue-matched trait. Asterisks indicate significance from 1-sided pairwise t-tests (*** P < 0.001, ** P < 0.01, * P 3) for each trait instead of the top 500 PoPS-prioritized disease gene set. Numerical results are provided in Supplementary Table S4 . Supplementary Figure S10. Additional results from the application of Spacelink on the mouse embryonic Stereo-seq SRT dataset . (a). Number of of SVGs identified by Spacelink at each stage within the non-brain region (top), and average ESV at each stage within the whole embryo (second), brain region (third), and non-brain region (bottom), with bars shaded from light to dark green across developmental time. Error bars denote 95% confidence intervals. (b). Volcano plots of gene-level regressions associating stage-corrected mean of sc-transform normalized gene expressions (adjusted for embryo size; Methods ) against developmental time in the whole embryo (left) and brain (right). (c). Volcano plots of gene-level regressions associating stage-corrected overall variance (spatial + non-spatial) of sc-transform normalized gene expressions (adjusted for embryo size) against developmental time in the whole embryo (left) and brain (right). In (b) and (c), each point represents the regression slope (x-axis) and –log10 p-value (y-axis). Genes with nominally significant (P < 0.05) negative and positive associations are shown in light and dark green, respectively. (d). Pathway enrichment analyses (based on ConsensusPathdb 63 ) of genes with significantly negative (light green) and positive (dark green) associations between mean (top panel) or variance (bottom panel) of sc-transform normalized gene expressions and developmental progression in the whole embryo (left column) and brain region (right column). Numerical results are reported in Supplementary Table S5 . Supplementary Figure S11. Odds ratios based on overlap of perturbation programs in each cell type for each pair of 35 autism knockout genes in in-vivo Perturb-seq assay . Odds ratios of excess overlap between downstream-altered genes (perturbation programs) in each of 5 cell types, and across cell types, for each pair of 35 de-novo autism knockout (KO) genes in Jin et al. 136 . The magnitude of overlap is represented by odds ratio (dot size), and significance by −log10(p-value) (dot color). Panels show results across all cell types (top left), astroglia (top middle), excitatory neurons (top right), inhibitory neurons (bottom left), microglia (bottom middle), and oligodendrocytes (ODCs; bottom right). Numerical results are reported in Supplementary Table S6 . Supplementary Figure S12. Additional results for Spacelink application on downstream perturbation programs in mouse in vivo perturbations of autism risk genes . (a). Effective Spatial Variability (ESV) of the 35 autism risk genes that were knocked out (KO) in the vivo Perturb-seq assay 136 across mouse organogenesis stages. (b). Average ESV of the KO genes (green) and downstream perturbation program genes in astroglia (dark blue), excitatory neurons (light blue), inhibitory neurons (red), microglia (yellow), and oligodendrocytes (orange) across organogenesis stages. Error bars denote 95% confidence intervals. (c). Average stage-corrected ESV of perturbation program genes in specific pairs of top three cell types with the strongest odds ratios of overlap in Figure 6b . Bars are shaded from light to dark green across developmental time. The left panel shows genes altered in astroglia and excitatory neurons but not microglia; the middle panel, in excitatory neurons and microglia but not astroglia; and the right panel, in astroglia and microglia but not excitatory neurons. Error bars denote 95% confidence intervals. (d) For each of the eight organogenesis stages, correlation between brain-specific and non–brain-specific mean (top) or variance (bottom) of sc-transform normalized gene expression across genes from three broad classes: (i) union of perturbation program genes in excitatory neurons and astroglia, (ii) union of perturbation program genes in inhibitory neurons, microglia, and oligodendrocytes, and (iii) all genes. (e). ESV of perturbation program genes in astroglia (top left), excitatory neurons (top right), inhibitory neurons (bottom left), microglia (bottom middle), and oligodendrocytes (bottom right) across eight organogenesis stages, stratified by the number of perturbed genes affecting them. Thick lines denote average ESV; shaded areas indicate confidence bands. The red dashed line shows the overall trend estimated by linear regression. Asterisks denote significance of Pearson correlation (*** P < 0.001). Numerical results are reported in Supplementary Table S6 . Supplementary Figure S13. Additional results integrating Spacelink with perturbation programs in mouse in vivo perturbations of autism risk genes . (a) Scatter plots of the average ESV of (i) all genes in genetically autism disease-associated gene set versus (ii) genes in the gene set that are in a perturbation program for some autism risk gene knockout (KO) in the in-vivo Perturb-seq assay. The disease-associated gene sets include (i) the top 100 MAGMA-prioritized polygenic genes (left), (ii) the top 100 PoPS-prioritized polygenic genes (middle), and (iii) 108 monogenic autism risk genes 73 , 74 (right). (b). Top: Scatter plots of the average mean of sc-transform normalized gene expression of (i) perturbation programs genes per perturbation in each stage versus (ii) subset of genes from (i) that overlap either of the three autism-associated gene sets. Bottom: Scatter plots of the overall mean of sc-transform normalized expression of (i) all genes in genetically autism disease-associated gene set versus (ii) genes in the gene set that are in a perturbation program for some autism risk gene knockout (KO) in the in-vivo Perturb-seq assay. (c). Top: Scatter plots of the average overall variance (spatial + non-spatial) of sc-transform normalized expression of (i) perturbation programs genes per perturbation in each stage versus (ii) subset of genes from (i) that overlap either of the three autism-associated gene sets. Bottom: Scatter plots of the average overall variance of sc-transform normalized expression of i) all genes in genetically autism disease-associated gene set versus (ii) genes in the gene set that are in a perturbation program for some autism risk gene knockout (KO) in the in-vivo Perturb-seq assay. (d). Density plots of perturbation program genes intersecting each autism-associated gene set, stratified by the number of perturbations affecting them. Shown are results across all cell types (top), as well as within astroglia (second), excitatory neurons (third), inhibitory neurons (fourth), microglia (fifth), and oligodendrocytes (bottom). Genes intersecting the top 100 MAGMA-prioritized polygenic genes are shown in red, those intersecting the top 100 PoPS-prioritized polygenic genes in blue, and those intersecting the 108 monogenic autism risk genes in yellow. In (a-c), red lines indicate the fitted regression line from a linear model with zero intercept; slopes are reported. Numerical results are reported in Supplementary Table S6 . Supplementary Figure S14. Additional results from the application of Spacelink to 10x Visium human DLPFC dataset from 32 samples at different stages of Alzheimer’s disease . (a). Boxplots of maximum normalized distances of tissue spots from each of 32 human samples in the 10x DLPFC dataset and 8 mouse samples in the Stereo-seq organogenesis dataset. For each sample, the maximum normalized distance was calculated by dividing all pairwise spot distances by the minimum distance and then taking the maximum of these normalized values. (b). Boxplots of the number of spatially variable genes (SVGs) (top panel) and the corrected or uncorrected average ESV (bottom panel) across genes at control and three stages of AD pathology, defined by amyloid burden (first column), tangle accumulation (second column), and neuritic plaque counts (third and fourth columns). For amyloid burden and tangle accumulation, CTRL includes 4 control samples from 2 healthy donors; AD-Lower includes 10 samples with the lowest values of the pathology indicator; AD-Mid includes 9 samples with intermediate values of the pathology indicator; and AD-Higher includes 9 samples with the highest values of the pathology indicator. For neuritic plaque counts, CTRL includes 4 control samples from 2 healthy donors; AD-Lower includes 8 samples with the lowest counts; AD-Mid includes 10 samples with intermediate counts; and AD-Higher includes 10 samples with the highest counts. The difference in sample numbers between the two classifications (4, 10, 9, 9 vs. 4, 8, 10, 10) arises because some samples had identical neuritic plaque counts, requiring their assignment to the same group. (c). Volcano plots of regression analyses associating ESV with neuritic plaque count at the gene level, adjusting for cortical layer composition, mean expression levels, and donor-level random effects. Each point represents the regression slope coefficient and its –log p-value for a gene. Genes nominally significant (P < 0.05) are highlighted in green. (d). Pathway enrichment analysis of 115 genes with nominally significant (P < 0.05) negative associations between spatial variability (ESV) and neuritic plaque count. Enrichment was performed using ConsensusPathDB 63 with the Visium panel genes as background. (e). Spatial expression plots of YWHAG (left column), with layer annotations (middle column) and plaque annotations (right column) in one control sample (top panel) and one high–tangle-accumulation sample (bottom panel). (f). Spatial expression plots of SNCB (left column), with layer annotations (middle column) and plaque annotations (right column) in one control sample (top panel) and one high–neuritic-plaque-count sample (bottom panel). (g). Volcano plots of regression analyses associating ESV with amyloid burden (left), tangle accumulation (middle), and neuritic plaque count (right) at the gene level, adjusting for cortical layer composition, mean expression levels, variance of expression levels, and donor-level random effects. Each point represents the regression slope coefficient and its –log p-value for a gene. Genes nominally significant (P < 0.05) are highlighted: amyloid (blue), tau burden (orange), and neuritic plaque counts (green). (h). Pathway enrichment analyses of 260, 166, and 103 genes with nominally significant (P < 0.05) negative associations between spatial variability (ESV) and amyloid burden (left), tangle accumulation (middle), and neuritic plaque count (right). Enrichment was performed using ConsensusPathDB with the Visium panel genes as background. Numerical results are reported in Supplementary Table S7 . Supplementary Figure S15. Additional results from the application of Spacelink to 10x Visium human DLPFC dataset from 32 samples at different stages of Alzheimer’s disease . (a). Volcano plots of regression analyses associating the overall mean of sc-transform normalized expression values with amyloid burden (left), tangle accumulation (middle), and neuritic plaque count (right) at the gene level, adjusting for cortical layer composition and donor-level random effects. Each point represents the regression slope coefficient and its –log p-value for a gene. Nominally significant genes (P < 0.05) are highlighted in blue (amyloid), orange (tau), and green (plaques). (b). Pathway enrichment analyses of 1,051 and 745 genes negatively associated with tangle accumulation (left) and neuritic plaque count (right) (P < 0.05), using ConsensusPathDB with Visium panel genes as background. (c). Boxplots of the number of SVGs (top) and average ESV (bottom) within each cortical layer in control and three AD groups defined by amyloid burden, tangle accumulation, or neuritic plaque counts. Group definitions (CTRL, AD-Lower, AD-Mid, AD-Higher) follow Supplementary Figure S15b . (d). Spatial plots of PGK1 (left) and GPI (right) in 10x Visium human DLPFC data (top) and mouse SRT data (bottom). For human data, expression, layer, and plaque annotations are shown for one control and one high–amyloid-burden sample. For mouse data, expression and layer annotations are shown for 4- and 12-month-old samples. (e). Spatial plots of MEF2C in human DLPFC (top) and mouse SRT (bottom). For human data, expression, layer, and plaque annotations are shown for one control and one high–amyloid-burden sample. For mouse data, expression and layer annotations are shown for 4- and 12-month-old samples. Numerical results are reported in Supplementary Table S7 . Acknowledgements K.K.D. was supported by R00HG012203 and R01HG014008, NIH/NCI Cancer Center Support Grant P30CA008748, AWS Imagine Fund and Josie Robertson Investigator Program. HUK and MT were supported by an Alzheimer’s Association Grant through the AD Strategic Fund (ADSF-21-816675). ROSMAP is supported by P30AG10161, P30AG72975, R01AG17917, R01AG015819, U01AG072572, and U01AG046152. Funder Information Declared National Human Genome Research Institute, https://ror.org/00baak391 , R00HG012203 , R01HG014008 National Cancer Institute, https://ror.org/040gcmg81 , P30CA008748 Amazon (United States), https://ror.org/04mv4n011 Memorial Sloan Kettering Cancer Center, https://ror.org/02yrq0923 Alzheimer’s Association, https://ror.org/0375f4d26 , ADSF-21-816675 National Institute of Neurological Disorders and Stroke, https://ror.org/01s5ya894 , P30AG10161 , P30AG72975 National Institute on Aging, https://ror.org/049v75w11 , R01AG17917 , R01AG015819 , U01AG072572 , U01AG046152 References 1. ↵ Ståhl , P. L. , Salmén , F. , Vickovic , S. , Lundmark , A. , Navarro , J. F. , Magnusson , J. , Giacomello , S. , Asp , M. , Westholm , J. O. , Huss , M. , Mollbrink , A. , Linnarsson , S. , Codeluppi , S. , Borg , Å. , Pontén , F. , Costea , P. I. , Sahlén , P. , Mulder , J. , Bergmann , O. , Lundeberg , J. & Frisén , J. Visualization and analysis of gene expression in tissue sections by spatial transcriptomics . Science 353 , 78 – 82 ( 2016 ). OpenUrl Abstract / FREE Full Text 2. Rao , A. , Barkley , D. , França , G. S. & Yanai , I. Exploring tissue architecture using spatial transcriptomics . Nature 596 , 211 – 220 ( 2021 ). OpenUrl CrossRef PubMed 3. ↵ Moses , L. & Pachter , L. Museum of spatial transcriptomics . Nat. Methods 19 , 534 – 546 ( 2022 ). OpenUrl CrossRef PubMed 4. ↵ Jain , S. & Eadon , M. T. Spatial transcriptomics in health and disease . Nat. Rev. Nephrol . 20 , 659 – 671 ( 2024 ). OpenUrl CrossRef PubMed 5. Williams , C. G. , Lee , H. J. , Asatsuma , T. , Vento-Tormo , R. & Haque , A. An introduction to spatial transcriptomics for biomedical research . Genome Med . 14 , 68 ( 2022 ). OpenUrl CrossRef PubMed 6. ↵ Miyoshi , E. , Morabito , S. , Henningfield , C. M. , Das , S. , Rahimzadeh , N. , Shabestari , S. K. , Michael , N. , Emerson , N. , Reese , F. , Shi , Z. , Cao , Z. , Srinivasan , S. S. , Scarfone , V. M. , Arreola , M. A. , Lu , J. , Wright , S. , Silva , J. , Leavy , K. , Lott , I. T. , Doran , E. , Yong , W. H. , Shahin , S. , Perez-Rosendahl , M. , Head , E. , Green , K. N. & Swarup , V. Spatial and single-nucleus transcriptomic analysis of genetic and sporadic forms of Alzheimer’s disease . Nat. Genet . 56 , 2704 – 2717 ( 2024 ). OpenUrl CrossRef PubMed 7. ↵ Song , L. , Chen , W. , Hou , J. , Guo , M. & Yang , J. Spatially resolved mapping of cells associated with human complex traits . Nature 641 , 932 – 941 ( 2025 ). OpenUrl CrossRef PubMed 8. ↵ Weber , L. M. , Saha , A. , Datta , A. , Hansen , K. D. & Hicks , S. C. nnSVG for the scalable identification of spatially variable genes using nearest-neighbor Gaussian processes . Nat. Commun . 14 , 4059 ( 2023 ). OpenUrl CrossRef PubMed 9. ↵ Svensson , V. , Teichmann , S. A. & Stegle , O. SpatialDE: identification of spatially variable genes . Nat. Methods 15 , 343 – 346 ( 2018 ). OpenUrl CrossRef PubMed 10. ↵ Kats , I. , Vento-Tormo , R. & Stegle , O. SpatialDE2: Fast and localized variance component analysis of spatial transcriptomics . 2021.10.27.466045 Preprint at doi: 10.1101/2021.10.27.466045 ( 2021 ) OpenUrl Abstract / FREE Full Text 11. ↵ Sun , S. , Zhu , J. & Zhou , X. Statistical analysis of spatial expression patterns for spatially resolved transcriptomic studies . Nat. Methods 17 , 193 – 200 ( 2020 ). OpenUrl CrossRef PubMed 12. ↵ Su , H. , Wu , Y. , Chen , B. & Cui , Y. STANCE: a unified statistical model to detect cell-type-specific spatially variable genes in spatial transcriptomics . Nat. Commun . 16 , 1793 ( 2025 ). OpenUrl PubMed 13. ↵ Shang , L. , Wu , P. & Zhou , X. Statistical identification of cell type-specific spatially variable genes in spatial transcriptomics . Nat. Commun . 16 , 1059 ( 2025 ). OpenUrl CrossRef PubMed 14. ↵ Chen , T.-Y. , You , L. , Hardillo , J. A. U. & Chien , M.-P. Spatial Transcriptomic Technologies . Cells 12 , 2042 ( 2023 ). OpenUrl CrossRef 15. ↵ Binan , L. , Jiang , A. , Danquah , S. A. , Valakh , V. , Simonton , B. , Bezney , J. , Manguso , R. T. , Yates , K. B. , Nehme , R. , Cleary , B. & Farhi , S. L. Simultaneous CRISPR screening and spatial transcriptomics reveal intracellular, intercellular, and functional transcriptional circuits . Cell 188 , 2141 – 2158.e18 ( 2025 ). OpenUrl CrossRef PubMed 16. ↵ Shen , K. , Seow , W. Y. , Keng , C. T. , Shern , D. L. , Guo , K. , Meliani , A. , Hajis , I. M. B. , Chen , K. H. & Chew , W. L. Spatial Perturb-Seq: Single-Cell Functional Genomics within Intact Tissue Architecture . 2024.12.19.628843 Preprint at doi: 10.1101/2024.12.19.628843 ( 2024 ) OpenUrl Abstract / FREE Full Text 17. ↵ Abedini , A. , Levinsohn , J. , Klötzer , K. A. , Dumoulin , B. , Ma , Z. , Frederick , J. , Dhillon , P. , Balzer , M. S. , Shrestha , R. , Liu , H. , Vitale , S. , Bergeson , A. M. , Devalaraja-Narashimha , K. , Grandi , P. , Bhattacharyya , T. , Hu , E. , Pullen , S. S. , Boustany-Kari , C. M. , Guarnieri , P. , Karihaloo , A. , Traum , D. , Yan , H. , Coleman , K. , Palmer , M. , Sarov-Blat , L. , Morton , L. , Hunter , C. A. , Kaestner , K. H. , Li , M. & Susztak , K. Single-cell multi-omic and spatial profiling of human kidneys implicates the fibrotic microenvironment in kidney disease progression . Nat. Genet . 56 , 1712 – 1724 ( 2024 ). OpenUrl CrossRef PubMed 18. Liu , Z. , Bian , X. , Luo , L. , Björklund , Å.K. , Li , L. , Zhang , L. , Chen , Y. , Guo , L. , Gao , J. , Cao , C. , Wang , J. , He , W. , Xiao , Y. , Zhu , L. , Annusver , K. , Gopee , N. H. , Basurto-Lozada , D. , Horsfall , D. , Bennett , C. L. , Kasper , M. , Haniffa , M. , Sommar , P. , Li , D. & Landén , N. X. Spatiotemporal single-cell roadmap of human skin wound healing . Cell Stem Cell 32 , 479 – 498.e8 ( 2025 ). OpenUrl CrossRef PubMed 19. ↵ Vannan , A. , Lyu , R. , Williams , A. L. , Negretti , N. M. , Mee , E. D. , Hirsh , J. , Hirsh , S. , Hadad , N. , Nichols , D. S. , Calvi , C. L. , Taylor , C. J. , Polosukhin , V. V. , Serezani , A. P. M. , McCall , A. S. , Gokey , J. J. , Shim , H. , Ware , L. B. , Bacchetta , M. J. , Shaver , C. M. , Blackwell , T. S. , Walia , R. , Sucre , J. M. S. , Kropski , J. A. , McCarthy , D. J. & Banovich , N. E. Spatial transcriptomics identifies molecular niche dysregulation associated with distal lung remodeling in pulmonary fibrosis . Nat. Genet . 57 , 647 – 658 ( 2025 ). OpenUrl CrossRef PubMed 20. ↵ Lawson , C. L. & Hanson , R. J. Solving Least Squares Problems . ( Society for Industrial and Applied Mathematics , 1995 ). doi: 10.1137/1.9781611971217 OpenUrl CrossRef 21. ↵ Cable , D. M. , Murray , E. , Zou , L. S. , Goeva , A. , Macosko , E. Z. , Chen , F. & Irizarry , R. A. Robust decomposition of cell type mixtures in spatial transcriptomics . Nat. Biotechnol . 40 , 517 – 526 ( 2022 ). OpenUrl CrossRef PubMed 22. ↵ Maestrini , L. , Hui , F. K. C. & Welsh , A. H. Restricted maximum likelihood estimation in generalized linear mixed models . Preprint at doi: 10.48550/arXiv.2402.12719 ( 2025 ) OpenUrl CrossRef 23. ↵ Pardo , B. , Spangler , A. , Weber , L. M. , Page , S. C. , Hicks , S. C. , Jaffe , A. E. , Martinowich , K. , Maynard , K. R. & Collado-Torres , L. spatialLIBD: an R/Bioconductor package to visualize spatially-resolved transcriptomics data . BMC Genomics 23 , 434 ( 2022 ). OpenUrl CrossRef PubMed 24. ↵ Andrews , T. S. , Nakib , D. , Perciani , C. T. , Ma , X. Z. , Liu , L. , Winter , E. , Camat , D. , Chung , S. W. , Lumanto , P. , Manuel , J. , Mangroo , S. , Hansen , B. , Arpinder , B. , Thoeni , C. , Sayed , B. , Feld , J. , Gehring , A. , Gulamhusein , A. , Hirschfield , G. M. , Ricciuto , A. , Bader , G. D. , McGilvray , I. D. & MacParland , S. Single-cell, single-nucleus, and spatial transcriptomics characterization of the immunological landscape in the healthy and PSC human liver . J. Hepatol . 80 , 730 – 743 ( 2024 ). OpenUrl CrossRef PubMed 25. ↵ Chen , A. , Liao , S. , Cheng , M. , Ma , K. , Wu , L. , Lai , Y. , Qiu , X. , Yang , J. , Xu , J. , Hao , S. , Wang , X. , Lu , H. , Chen , X. , Liu , X. , Huang , X. , Li , Z. , Hong , Y. , Jiang , Y. , Peng , J. , Liu , S. , Shen , M. , Liu , C. , Li , Q. , Yuan , Y. , Wei , X. , Zheng , H. , Feng , W. , Wang , Z. , Liu , Y. , Wang , Z. , Yang , Y. , Xiang , H. , Han , L. , Qin , B. , Guo , P. , Lai , G. , Muñoz-Cánoves , P. , Maxwell , P. H. , Thiery , J. P. , Wu , Q.-F. , Zhao , F. , Chen , B. , Li , M. , Dai , X. , Wang , S. , Kuang , H. , Hui , J. , Wang , L. , Fei , J.-F. , Wang , O. , Wei , X. , Lu , H. , Wang , B. , Liu , S. , Gu , Y. , Ni , M. , Zhang , W. , Mu , F. , Yin , Y. , Yang , H. , Lisby , M. , Cornall , R. J. , Mulder , J. , Uhlén , M. , Esteban , M. A. , Li , Y. , Liu , L. , Xu , X. & Wang , J. Spatiotemporal transcriptomic atlas of mouse organogenesis using DNA nanoball-patterned arrays . Cell 185 , 1777 – 1792.e21 ( 2022 ). OpenUrl CrossRef PubMed 26. ↵ Karaahmet , B. , Zhang , Y. , Duquesne , L. , Sigalov , A. , Yung , C. , Kroshilina , A. , Bennett , D. A. , Taga , M. & Klein , H.-U. Spatial Transcriptomic Analysis Identifies a SERPINA3 - Expressing Astrocytic State Associated with the Human Neuritic Plaque Microenvironment . Preprint at doi: 10.1101/2024.11.13.623438 ( 2024 ) OpenUrl Abstract / FREE Full Text 27. ↵ Weeks , E. M. , Ulirsch , J. C. , Cheng , N. Y. , Trippe , B. L. , Fine , R. S. , Miao , J. , Patwardhan , T. A. , Kanai , M. , Nasser , J. , Fulco , C. P. , Tashman , K. C. , Aguet , F. , Li , T. , Ordovas-Montanes , J. , Smillie , C. S. , Biton , M. , Shalek , A. K. , Ananthakrishnan , A. N. , Xavier , R. J. , Regev , A. , Gupta , R. M. , Lage , K. , Ardlie , K. G. , Hirschhorn , J. N. , Lander , E. S. , Engreitz , J. M. & Finucane , H. K. Leveraging polygenic enrichments of gene features to predict genes underlying complex traits and diseases . Nat. Genet . 55 , 1267 – 1276 ( 2023 ). OpenUrl CrossRef PubMed 28. ↵ Leeuw , C. A. de , Mooij , J. M. , Heskes , T. & Posthuma , D. MAGMA: Generalized Gene-Set Analysis of GWAS Data . PLOS Comput. Biol . 11 , e1004219 ( 2015 ). OpenUrl CrossRef PubMed 29. ↵ de Leeuw , C. A. , Neale , B. M. , Heskes , T. & Posthuma , D. The statistical properties of gene-set analysis . Nat. Rev. Genet . 17 , 353 – 364 ( 2016 ). OpenUrl CrossRef PubMed 30. ↵ Moran , P. a. P. Notes on continuous stochastic phenomena . Biometrika 37 , 17 – 23 ( 1950 ). OpenUrl CrossRef PubMed Web of Science 31. ↵ SPARK-X: non-parametric modeling enables scalable and robust detection of spatial expression patterns for large spatial transcriptomic studies | Genome Biology | Full Text . At 32. ↵ Li , Z.M. , Patel , Z. , Song , D. , Yasa , S. N. , Cannoodt , R. , Yan , G. , Li , J. J. & Pinello , L. Systematic benchmarking of computational methods to identify spatially variable genes . Genome Biol . 26 , 285 ( 2025 ). OpenUrl PubMed 33. ↵ Song , D. , Wang , Q. , Yan , G. , Liu , T. , Sun , T. & Li , J. J. scDesign3 generates realistic in silico data for multimodal single-cell and spatial omics . Nat. Biotechnol . 42 , 247 – 252 ( 2024 ). OpenUrl CrossRef PubMed 34. ↵ Chen , C. , Kim , H. J. & Yang , P. Evaluating spatially variable gene detection methods for spatial transcriptomics data . Genome Biol . 25 , 18 ( 2024 ). OpenUrl CrossRef PubMed 35. ↵ Dey , K. K. , Hsiao , C. J. & Stephens , M. Visualizing the structure of RNA-seq expression data using grade of membership models . PLOS Genet . 13 , e1006599 ( 2017 ). OpenUrl CrossRef PubMed 36. ↵ Datta , A. , Banerjee , S. , Finley , A. O. & Gelfand , A. E. Hierarchical Nearest-Neighbor Gaussian Process Models for Large Geostatistical Datasets . J. Am. Stat. Assoc . 111 , 800 – 812 ( 2016 ). OpenUrl CrossRef PubMed 37. ↵ Pardo , B. , Spangler , A. , Weber , L. M. , Page , S. C. , Hicks , S. C. , Jaffe , A. E. , Martinowich , K. , Maynard , K. R. & Collado-Torres , L. spatialLIBD: an R/Bioconductor package to visualize spatially-resolved transcriptomics data . BMC Genomics 23 , 434 ( 2022 ). OpenUrl CrossRef PubMed 38. ↵ Wang , X. , Allen , W. E. , Wright , M. A. , Sylwestrak , E. L. , Samusik , N. , Vesuna , S. , Evans , K. , Liu , C. , Ramakrishnan , C. , Liu , J. , Nolan , G. P. Bava , F.-A. & Deisseroth , K. Three-dimensional intact-tissue sequencing of single-cell transcriptional states . Science 361 , eaat5691 ( 2018 ). OpenUrl Abstract / FREE Full Text 39. ↵ Jagadeesh , K. A. , Dey , K. K. , Montoro , D. T. , Mohan , R. , Gazal , S. , Engreitz , J. M. , Xavier , R. J. , Price , A.L. & Regev , A. Identifying disease-critical cell types and cellular processes by integrating single-cell RNA-sequencing and human genetics . Nat. Genet . 54 , 1479 – 1492 ( 2022 ). OpenUrl CrossRef PubMed 40. Zhang , M. J. , Hou , K. , Dey , K. K. , Sakaue , S. , Jagadeesh , K. A. , Weinand , K. , Taychameekiatchai , A. , Rao , P. , Pisco , A. O. , Zou , J. , Wang , B. , Gandal , M. , Raychaudhuri , S. , Pasaniuc , B. & Price , A. L. Polygenic enrichment distinguishes disease associations of individual cells in single-cell RNA-seq data . Nat. Genet . 54 , 1572 – 1580 ( 2022 ). OpenUrl CrossRef PubMed 41. ↵ Fabiha , T. , Raine , I. , Kundu , S. , Pampari , A. , Abramov , S. , Boytsov , A. , Strouse , K. , Dura , K. , Fang , W. , Kerner , G. , Butts , J. , Ali , T. , Gschwind , A. , Mualim , K. S. , Moore , J. E. , Weng , Z. , Ulirsch , J. , Ji , H. E. , Vierstra , J. , Reddy , T. E. , Montgomery , S. B. , Engreitz , J. , Kundaje , A. , Tewhey , R. , Price , A. & Dey , K. A consensus variant-to-function score to functionally prioritize variants for disease . 2024.11.07.622307 Preprint at doi: 10.1101/2024.11.07.622307 ( 2024 ) OpenUrl Abstract / FREE Full Text 42. ↵ Aihara , G. , Clifton , K. , Chen , M. , Li , Z. , Atta , L. , Miller , B. F. , Satija , R. , Hickey , J. W. & Fan , J. SEraster: a rasterization preprocessing framework for scalable spatial omics data analysis . Bioinformatics 40 , btae412 ( 2024 ). OpenUrl PubMed 43. ↵ Jagadeesh , K. A. , Dey , K. K. , Montoro , D. T. , Mohan , R. , Gazal , S. , Engreitz , J. M. , Xavier , R. J. , Price , A. L. & Regev , A. Identifying disease-critical cell types and cellular processes by integrating single-cell RNA sequencing and human genetics . Nat. Genet . 54 , 1479 – 1492 ( 2022 ). OpenUrl CrossRef PubMed 44. ↵ Finucane , H. K. , Bulik-Sullivan , B. , Gusev , A. , Trynka , G. , Reshef , Y. , Loh , P.-R. , Anttila , V. , Xu , H. , Zang , C. , Farh , K. , Ripke , S. , Day , F. R. , ReproGen Consortium , Schizophrenia Working Group of the Psychiatric Genomics Consortium , RACI Consortium , Purcell , S. , Stahl , E. , Lindstrom , S. , Perry , J. R. B. , Okada , Y. , Raychaudhuri , S. , Daly , M. J. , Patterson , N. , Neale , B. M. & Price , A. L. Partitioning heritability by functional annotation using genome-wide association summary statistics . Nat. Genet . 47 , 1228 – 1235 ( 2015 ). OpenUrl CrossRef PubMed 45. ↵ Jagadeesh , K. A. , Dey , K. K. , Montoro , D. T. , Mohan , R. , Gazal , S. , Engreitz , J. M. , Xavier , R. J. , Price , A. L. & Regev , A. Identifying disease-critical cell types and cellular processes by integrating single-cell RNA-sequencing and human genetics . Nat. Genet . 54 , 1479 – 1492 ( 2022 ). OpenUrl CrossRef PubMed 46. ↵ Howson , J. M. M. , Zhao , W. , Barnes , D. R. , Ho , W.-K. , Young , R. , Paul , D. S. , Waite , L. L. , Freitag , D. F. , Fauman , E. B. , Salfati , E. L. , Sun , B. B. , Eicher , J. D. , Johnson , A. D. , Sheu , W. H. H. , Nielsen , S. F. , Lin , W.-Y. , Surendran , P. , Malarstig , A. , Wilk , J. B. , Tybjærg-Hansen , A. , Rasmussen , K. L. , Kamstrup , P. R. , Deloukas , P. , Erdmann , J. , Kathiresan , S. , Samani , N. J. , Schunkert , H. , Watkins , H. , CARDIoGRAMplusC4D , Do , R. , Rader , D. J. , Johnson , J. A. , Hazen , S. L. , Quyyumi , A. A. , Spertus , J. A. , Pepine , C. J. , Franceschini , N. , Justice , A. , Reiner , A. P. , Buyske , S. , Hindorff , L. A. , Carty , C. L. , North , K. E. , Kooperberg , C. , Boerwinkle , E. , Young , K. , Graff , M. , Peters , U. , Absher , D. , Hsiung , C. A. , Lee , W.-J. , Taylor , K. D. , Chen , Y.-H. , Lee , I.-T. , Guo , X. , Chung , R.-H. , Hung , Y.-J. , Rotter , J. I. , Juang , J.-M. J. , Quertermous , T. , Wang , T.-D. , Rasheed , A. , Frossard , P. , Alam , D. S. , Majumder , A. A. S. , Di Angelantonio , E. , Chowdhury , R. , EPIC-CVD , Chen , Y.-D.I. , Nordestgaard , B. G. , Assimes , T. L. , Danesh , J. , Butterworth , A. S. & Saleheen , D. Fifteen new risk loci for coronary artery disease highlight arterial-wall-specific mechanisms . Nat. Genet . 49 , 1113 – 1119 ( 2017 ). OpenUrl CrossRef PubMed 47. ↵ Bycroft , C. , Freeman , C. , Petkova , D. , Band , G. , Elliott , L. T. , Sharp , K. , Motyer , A. , Vukcevic , D. , Delaneau , O. , O’Connell , J. , Cortes , A. , Welsh , S. , Young , A. , Effingham , M. , McVean , G. , Leslie , S. , Allen , N. , Donnelly , P. & Marchini , J. The UK Biobank resource with deep phenotyping and genomic data . Nature 562 , 203 – 209 ( 2018 ). OpenUrl CrossRef PubMed 48. ↵ Bipolar Disorder and Schizophrenia Working Group of the Psychiatric Genomics Consortium . Electronic address: douglas.ruderfer{at}vanderbilt.edu & Bipolar Disorder and Schizophrenia Working Group of the Psychiatric Genomics Consortium. Genomic Dissection of Bipolar Disorder and Schizophrenia, Including 28 Subphenotypes . Cell 173 , 1705 – 1715.e16 ( 2018 ). OpenUrl CrossRef PubMed 49. ↵ Pepin , M. E. & Gupta , R. M. The Role of Endothelial Cells in Atherosclerosis: Insights from Genetic Association Studies . Am. J. Pathol . 194 , 499 – 509 ( 2024 ). OpenUrl CrossRef PubMed 50. ↵ Matsuzawa , Y. & Lerman , A. Endothelial Dysfunction and Coronary Artery Disease: Assessment, Prognosis and Treatment . Coron. Artery Dis . 25 , 713 – 724 ( 2014 ). OpenUrl CrossRef PubMed 51. ↵ Hsu , J.-H. , Chien , I.-C. , Lin , C.-H. Chou , Y.-J. & Chou , P. Increased risk of chronic liver disease in patients with schizophrenia: a population-based cohort study . Psychosomatics 55 , 163 – 171 ( 2014 ). OpenUrl PubMed 52. ↵ Grant , R. K. , Brindle , W. M. , Donnelly , M. C. , McConville , P. M. , Stroud , T. G. , Bandieri , L. & Plevris , J. N. Gastrointestinal and liver disease in patients with schizophrenia: A narrative review . World J. Gastroenterol . 28 , 5515 – 5529 ( 2022 ). OpenUrl PubMed 53. ↵ Rogalski , J. & Tomczak , T. Ferroptosis as a potential connection between schizophrenia and metabolic dysfunction-associated steatotic liver disease – a narrative review . Adv. Psychiatry Neurol . 33 , 178 – 187 ( 2024 ). OpenUrl 54. ↵ Vannan , A. , Lyu , R. , Williams , A. L. , Negretti , N. M. , Mee , E. D. , Hirsh , J. , Hirsh , S. , Hadad , N. , Nichols , D. S. , Calvi , C. L. , Taylor , C. J. , Polosukhin , V. V. , Serezani , A. P. M. , McCall , A. S. , Gokey , J. J. , Shim , H. , Ware , L. B. , Bacchetta , M. J. , Shaver , C. M. , Blackwell , T. S. , Walia , R. , Sucre , J. M. S. , Kropski , J. A. , McCarthy , D. J. & Banovich , N. E. Spatial transcriptomics identifies molecular niche dysregulation associated with distal lung remodeling in pulmonary fibrosis . Nat. Genet . 57 , 647 – 658 ( 2025 ). OpenUrl CrossRef PubMed 55. ↵ Long , H. , Lichtnekert , J. , Andrassy , J. , Schraml , B. U. , Romagnani , P. & Anders , H.-J. Macrophages and fibrosis: how resident and infiltrating mononuclear phagocytes account for organ injury, regeneration or atrophy . Front. Immunol . 14 , ( 2023 ). 56. ↵ Mathys , H. , Davila-Velderrain , J. , Peng , Z. , Gao , F. , Mohammadi , S. , Young , J. Z. , Menon , M. , He , L. , Abdurrob , F. , Jiang , X. , Martorell , A. J. , Ransohoff , R. M. , Hafler , B. P. , Bennett , D. A. , Kellis , M. & Tsai , L.-H. Single-cell transcriptomic analysis of Alzheimer’s disease . Nature 570 , 332 – 337 ( 2019 ). OpenUrl CrossRef PubMed 57. ↵ Lee , S.-E. & Lee , G. H. Reelin Affects Signaling Pathways of a Group of Inhibitory Neurons and the Development of Inhibitory Synapses in Primary Neurons . Int. J. Mol. Sci . 22 , 7510 ( 2021 ). OpenUrl PubMed 58. ↵ Okada , Y. , Wu , D. , Trynka , G. , Raj , T. , Terao , C. , Ikari , K. , Kochi , Y. , Ohmura , K. , Suzuki , A. , Yoshida , S. , Graham , R. R. , Manoharan , A. , Ortmann , W. , Bhangale , T. , Denny , J. C. , Carroll , R. J. , Eyler , A. E. , Greenberg , J. D. , Kremer , J. M. , Pappas , D. A. , Jiang , L. , Yin , J. , Ye , L. , Su , D.-F. , Yang , J. , Xie , G. , Keystone , E. , Westra , H.-J. , Esko , T. , Metspalu , A. , Zhou , X. , Gupta , N. , Mirel , D. , Stahl , E. A. , Diogo , D. , Cui , J. , Liao , K. , Guo , M. H. , Myouzen , K. , Kawaguchi , T. , Coenen , M. J. H. , van Riel , P. L. C. M. , van de Laar , M. A. F. J. , Guchelaar , H.-J. , Huizinga , T. W. J. , Dieudé , P. , Mariette , X. , Louis Bridges Jr , S. , Zhernakova , A. , Toes , R. E. M. , Tak , P. P. , Miceli-Richard , C. , Bang , S.-Y. , Lee , H.-S. , Martin , J. , Gonzalez-Gay , M. A. , Rodriguez-Rodriguez , L. , Rantapää-Dahlqvist , S. , Ärlestig , L. , Choi , H. K. , Kamatani , Y. , Galan , P. , Lathrop , M. , Eyre , S. , Bowes , J. , Barton , A. , de Vries , N. , Moreland , L. W. , Criswell , L. A. , Karlson , E. W. , Taniguchi , A. , Yamada , R. , Kubo , M. , Liu , J. S. , Bae , S.-C. , Worthington , J. , Padyukov , L. , Klareskog , L. , Gregersen , P. K. , Raychaudhuri , S. , Stranger , B. E. , De Jager , P. L. , Franke , L. , Visscher , P. M. , Brown , M. A. , Yamanaka , H. , Mimori , T. , Takahashi , A. , Xu , H. , Behrens , T. W. , Siminovitch , K. A. , Momohara , S. , Matsuda , F. , Yamamoto , K. & Plenge , R. M. Genetics of rheumatoid arthritis contributes to biology and drug discovery . Nature 506 , 376 – 381 ( 2014 ). OpenUrl CrossRef PubMed Web of Science 59. ↵ Pappworth , I. Y. , Kulik , L. , Haluszczak , C. , Reuter , J. W. , Holers , V. M. & Marchbank , K. J. Increased B cell deletion and significantly reduced auto-antibody titre due to premature expression of human complement receptor 2 (CR2, CD21) . Mol. Immunol . 46 , 1042 – 1049 ( 2009 ). OpenUrl CrossRef PubMed 60. ↵ Mosca , A. , Braghini , M. R. , Andolina , G. , De Stefanis , C. , Cesarini , L. , Pastore , A. , Comparcola , D. , Monti , L. , Francalanci , P. , Balsano , C. , Pietrobattista , A. , Alisi , A. & Panera , N. Levels of Growth Differentiation Factor 15 Correlated with Metabolic Dysfunction-Associated Steatotic Liver Disease in Children . Int. J. Mol. Sci . 26 , 6486 ( 2025 ). OpenUrl PubMed 61. ↵ Galuppo , B. , Agazzi , C. , Pierpont , B. , Chick , J. , Li , Z. , Caprio , S. & Santoro , N. Growth differentiation factor 15 (GDF15) is associated with non-alcoholic fatty liver disease (NAFLD) in youth with overweight or obesity . Nutr. Diabetes 12 , 9 ( 2022 ). OpenUrl CrossRef PubMed 62. ↵ Hafemeister , C. & Satija , R. Normalization and variance stabilization of single-cell RNA-seq data using regularized negative binomial regression . Genome Biol . 20 , 296 ( 2019 ). OpenUrl CrossRef PubMed 63. ↵ Kamburov , A. & Herwig , R. ConsensusPathDB 2022: molecular interactions update as a resource for network biology . Nucleic Acids Res . 50 , D587 – D595 ( 2021 ). OpenUrl 64. ↵ Takahashi , H. & Liu , F.-C. Genetic patterning of the mammalian telencephalon by morphogenetic molecules and transcription factors . Birth Defects Res. Part C Embryo Today Rev . 78 , 256 – 266 ( 2006 ). OpenUrl 65. Green , D. G. , Whitener , A. E. , Mohanty , S. , Mistretta , B. , Gunaratne , P. , Yeh , A. T. & Lekven , A. C. Wnt signaling regulates neural plate patterning in distinct temporal phases with dynamic transcriptional outputs . Dev. Biol . 462 , 152 – 164 ( 2020 ). OpenUrl PubMed 66. ↵ Wu , J. , Yang , J. & Klein , P. S. Neural crest induction by the canonical Wnt pathway can be dissociated from anterior-posterior neural patterning in Xenopus . Dev. Biol . 279 , 220 – 232 ( 2005 ). OpenUrl CrossRef PubMed 67. ↵ Wang , Y. , Dye , C. A. , Sohal , V. , Long , J. E. , Estrada , R. C. , Roztocil , T. , Lufkin , T. , Deisseroth , K. , Baraban , S. C. & Rubenstein , J. L. R. Dlx5 and Dlx6 Regulate the Development of Parvalbumin-Expressing Cortical Interneurons . J. Neurosci . 30 , 5334 – 5345 ( 2010 ). OpenUrl Abstract / FREE Full Text 68. ↵ Robledo , R. F. , Rajan , L. , Li , X. & Lufkin , T. The Dlx5 and Dlx6 homeobox genes are essential for craniofacial, axial, and appendicular skeletal development . Genes Dev . 16 , 1089 – 1101 ( 2002 ). OpenUrl Abstract / FREE Full Text 69. ↵ Andersson , L. S. , Larhammar , M. , Memic , F. , Wootz , H. , Schwochow , D. , Rubin , C.-J. , Patra , K. , Arnason , T. , Wellbring , L. , Hjälm , G. , Imsland , F. , Petersen , J. L. , McCue , M. E. , Mickelson , J. R. , Cothran , G. , Ahituv , N. , Roepstorff , L. , Mikko , S. , Vallstedt , A. , Lindgren , G. , Andersson , L. & Kullander , K. Mutations in DMRT3 affect locomotion in horses and spinal circuit function in mice . Nature 488 , 642 – 646 ( 2012 ). OpenUrl CrossRef PubMed Web of Science 70. ↵ Perry , S. , Larhammar , M. , Vieillard , J. , Nagaraja , C. , Hilscher , M. M. , Tafreshiha , A. , Rofo , F. , Caixeta , F. V. & Kullander , K. Characterization of Dmrt3-Derived Neurons Suggest a Role within Locomotor Circuits . J. Neurosci . 39 , 1771 – 1782 ( 2019 ). OpenUrl Abstract / FREE Full Text 71. ↵ Lewis , M. T. Homeobox genes in mammary gland development and neoplasia . Breast Cancer Res . 2 , 158 ( 2000 ). OpenUrl CrossRef PubMed Web of Science 72. ↵ Taniguchi , Y. , Tanaka , O. , Sekiguchi , M. , Takekoshi , S. , Tsukamoto , H. , Kimura , M. , Imai , K. & Inoko , H. Enforced expression of the transcription factor HOXD3 under the control of the Wnt1 regulatory element modulates cell adhesion properties in the developing mouse neural tube . J. Anat . 219 , 589 – 600 ( 2011 ). OpenUrl PubMed 73. ↵ Amberger , J. S. , Bocchini , C. A. , Schiettecatte , F. , Scott , A. F. & Hamosh , A. OMIM.org: Online Mendelian Inheritance in Man (OMIM®), an online catalog of human genes and genetic disorders . Nucleic Acids Res . 43 , D789 – 798 ( 2015 ). OpenUrl CrossRef PubMed 74. ↵ Freund , M. K. , Burch , K. S. , Shi , H. , Mancuso , N. , Kichaev , G. , Garske , K. M. , Pan , D. Z. , Miao , Z. , Mohlke , K. L. , Laakso , M. , Pajukanta , P. , Pasaniuc , B. & Arboleda , V. A. Phenotype-Specific Enrichment of Mendelian Disorder Genes near GWAS Regions across 62 Complex Traits . Am. J. Hum. Genet . 103 , 535 – 552 ( 2018 ). OpenUrl CrossRef PubMed 75. ↵ Jin , X. , Simmons , S. K. , Guo , A. , Shetty , A. S. , Ko , M. , Nguyen , L. , Jokhi , V. , Robinson , E. , Oyler , P. , Curry , N. , Deangeli , G. , Lodato , S. , Levin , J. Z. , Regev , A. , Zhang , F. & Arlotta , P. In vivo Perturb-Seq reveals neuronal and glial abnormalities associated with autism risk genes . Science 370 , eaaz6063 ( 2020 ). OpenUrl Abstract / FREE Full Text 76. ↵ Satterstrom , F. K. , Kosmicki , J. A. , Wang , J. , Breen , M. S. , De Rubeis , S. , An , J.-Y. , Peng , M. , Collins , R. , Grove , J. , Klei , L. , Stevens , C. , Reichert , J. , Mulhern , M. S. , Artomov , M. , Gerges , S. , Sheppard , B. , Xu , X. , Bhaduri , A. , Norman , U. , Brand , H. , Schwartz , G. , Nguyen , R. , Guerrero , E. E. , Dias , C. , Aleksic , B. , Anney , R. , Barbosa , M. , Bishop , S. , Brusco , A. , Bybjerg-Grauholm , J. , Carracedo , A. , Chan , M. C. Y. , Chiocchetti , A. G. , Chung , B. H. Y. , Coon , H. , Cuccaro , M. L. , Curró , A. , Dalla Bernardina , B. , Doan , R. , Domenici , E. , Dong , S. , Fallerini , C. , Fernández-Prieto , M. , Ferrero , G. B. , Freitag , C. M. , Fromer , M. , Gargus , J. J. , Geschwind , D. , Giorgio , E. , González-Peñas , J. , Guter , S. , Halpern , D. , Hansen-Kiss , E. , He , X. , Herman , G. E. , Hertz-Picciotto , I. , Hougaard , D. M. , Hultman , C. M. , Ionita-Laza , I. , Jacob , S. , Jamison , J. , Jugessur , A. , Kaartinen , M. , Knudsen , G. P. , Kolevzon , A. , Kushima , I. , Lee , S. L. , Lehtimäki , T. , Lim , E. T. , Lintas , C. , Lipkin , W. I. , Lopergolo , D. , Lopes , F. , Ludena , Y. , Maciel , P. , Magnus , P. , Mahjani , B. , Maltman , N. , Manoach , D. S. , Meiri , G. , Menashe , I. , Miller , J. , Minshew , N. , Montenegro , E. M. S. , Moreira , D. , Morrow , E. M. , Mors , O. , Mortensen , P. B. , Mosconi , M. , Muglia , P. , Neale , B. M. , Nordentoft , M. , Ozaki , N. , Palotie , A. , Parellada , M. , Passos-Bueno , M. R. , Pericak-Vance , M. , Persico , A. M. , Pessah , I. , Puura , K. , Reichenberg , A. , Renieri , A. , Riberi , E. , Robinson , E. B. , Samocha , K. E. , Sandin , S. , Santangelo , S. L. , Schellenberg , G. , Scherer , S. W. , Schlitt , S. , Schmidt , R. , Schmitt , L. , Silva , I. M. W. , Singh , T. , Siper , P. M. , Smith , M. , Soares , G. , Stoltenberg , C. , Suren , P. , Susser , E. , Sweeney , J. , Szatmari , P. , Tang , L. , Tassone , F. , Teufel , K. , Trabetti , E. , Trelles , M. del P. , Walsh , C. A. , Weiss , L. A. , Werge , T. , Werling , D. M. , Wigdor , E. M. , Wilkinson , E. , Willsey , A. J. , Yu , T. W. , Yu , M. H. C. , Yuen , R. , Zachi , E. , Agerbo , E. , Als , T. D. , Appadurai , V. , Bækvad-Hansen , M. , Belliveau , R. , Buil , A. , Carey , C. E. , Cerrato , F. , Chambert , K. , Churchhouse , C. , Dalsgaard , S. , Demontis , D. , Dumont , A. , Goldstein , J. , Hansen , C. S. , Hauberg , M. E. , Hollegaard , M. V. , Howrigan , D. P. , Huang , H. , Maller , J. , Martin , A. R. , Martin , J. , Mattheisen , M. , Moran , J. , Pallesen , J. , Palmer , D. S. , Pedersen , C. B. , Pedersen , M. G. , Poterba , T. , Poulsen , J. B. , Ripke , S. , Schork , A. J. , Thompson , W. K. , Turley , P. , Walters , R. K. , Betancur , C. , Cook , E. H. , Gallagher , L. , Gill , M. , Sutcliffe , J. S. , Thurm , A. , Zwick , M. E. , Børglum , A. D. , State , M. W. , Cicek , A. E. , Talkowski , M. E. , Cutler , D. J. , Devlin , B. , Sanders , S. J. , Roeder , K. , Daly , M. J. & Buxbaum , J. D. Large-Scale Exome Sequencing Study Implicates Both Developmental and Functional Changes in the Neurobiology of Autism . Cell 180 , 568 – 584.e23 ( 2020 ). OpenUrl CrossRef PubMed 77. ↵ Sanders , S. J. , Murtha , M. T. , Gupta , A. R. , Murdoch , J. D. , Raubeson , M. J. , Willsey , A. J. , Ercan-Sencicek , A. G. , DiLullo , N. M. , Parikshak , N. N. , Stein , J. L. , Walker , M. F. , Ober , G. T. , Teran , N. A. , Song , Y. , El-Fishawy , P. , Murtha , R. C. , Choi , M. , Overton , J. D. , Bjornson , R. D. , Carriero , N. J. , Meyer , K. A. , Bilguvar , K. , Mane , S. M. , Šestan , N. , Lifton , R. P. , Günel , M. , Roeder , K. , Geschwind , D. H. , Devlin , B. & State , M. W. De novo mutations revealed by whole exome sequencing are strongly associated with autism . Nature 485 , 237 – 241 ( 2012 ). OpenUrl CrossRef PubMed Web of Science 78. ↵ Geiger-Schuller , K. , Eraslan , B. , Kuksenko , O. , Dey , K. K. , Jagadeesh , K. A. , Thakore , P. I. , Karayel , O. , Yung , A. R. , Rajagopalan , A. , Meireles , A. M. , Yang , K. D. , Amir-Zilberstein , L. , Delorey , T. , Phillips , D. , Raychowdhury , R. , Moussion , C. , Price , A. L. , Hacohen , N. , Doench , J. G. , Uhler , C. , Rozenblatt-Rosen , O. & Regev , A. Systematically characterizing the roles of E3-ligase family members in inflammatory responses with massively parallel Perturb-seq . 2023.01.23.525198 Preprint at doi: 10.1101/2023.01.23.525198 ( 2023 ) OpenUrl Abstract / FREE Full Text 79. ↵ Coquand , L. , Brunet Avalos , C. , Macé , A.-S. , Farcy , S. , Di Cicco , A. , Lampic , M. , Wimmer , R. , Bessières , B. , Attie-Bitach , T. , Fraisier , V. , Sens , P. , Guimiot , F. Brault , J.-B. & Baffet , A. D. A cell fate decision map reveals abundant direct neurogenesis bypassing intermediate progenitors in the human developing neocortex . Nat. Cell Biol . 26 , 698 – 709 ( 2024 ). OpenUrl CrossRef PubMed 80. ↵ Miranda-Negrón , Y. & García-Arrarás , J. E. Radial glia and radial glia-like cells: Their role in neurogenesis and regeneration . Front. Neurosci . 16 , 1006037 ( 2022 ). OpenUrl CrossRef PubMed 81. ↵ Allen , N. J. & Eroglu , C. Cell biology of astrocyte-synapse interactions . Neuron 96 , 697 – 708 ( 2017 ). OpenUrl CrossRef PubMed 82. ↵ Ricci , G. , Volpi , L. , Pasquali , L. , Petrozzi , L. & Siciliano , G. Astrocyte–neuron interactions in neurological disorders . J. Biol. Phys . 35 , 317 – 336 ( 2009 ). OpenUrl CrossRef PubMed 83. ↵ Replogle , J. M. , Saunders , R. A. , Pogson , A. N. , Hussmann , J. A. , Lenail , A. , Guna , A. , Mascibroda , L. , Wagner , E. J. , Adelman , K. , Lithwick-Yanai , G. , Iremadze , N. , Oberstrass , F. , Lipson , D. , Bonnar , J. L. , Jost , M. , Norman , T. M. & Weissman , J. S. Mapping information-rich genotype-phenotype landscapes with genome-scale Perturb-seq . Cell 185 , 2559 – 2575.e28 ( 2022 ). OpenUrl CrossRef PubMed 84. ↵ Joung , J. , Ma , S. , Tay , T. , Geiger-Schuller , K. R. , Kirchgatterer , P. C. , Verdine , V. K. , Guo , B. , Arias-Garcia , M. A. , Allen , W. E. , Singh , A. , Kuksenko , O. , Abudayyeh , O. O. , Gootenberg , J. S. , Fu , Z. , Macrae , R. K. , Buenrostro , J. D. , Regev , A. & Zhang , F. A transcription factor atlas of directed differentiation . Cell 186 , 209 – 229.e26 ( 2023 ). OpenUrl CrossRef PubMed 85. ↵ Matcovitch-Natan , O. , Winter , D. R. , Giladi , A. , Vargas Aguilar , S. , Spinrad , A. , Sarrazin , S. , Ben-Yehuda , H. , David , E. , Zelada González , F. , Perrin , P. , Keren-Shaul , H. , Gury , M. , Lara-Astaiso , D. , Thaiss , C. A. , Cohen , M. , Bahar Halpern , K. , Baruch , K. , Deczkowska , A. , Lorenzo-Vivas , E. , Itzkovitz , S. , Elinav , E. , Sieweke , M. H. , Schwartz , M. & Amit , I. Microglia development follows a stepwise program to regulate brain homeostasis . Science 353 , aad8670 ( 2016 ). OpenUrl Abstract / FREE Full Text 86. ↵ Winkler , C. & Yao , S. The midkine family of growth factors: diverse roles in nervous system formation and maintenance . Br. J. Pharmacol . 171 , 905 – 912 ( 2014 ). OpenUrl CrossRef PubMed 87. ↵ Campbell , W. A. , Fritsch-Kelleher , A. , Palazzo , I. , Hoang , T. , Blackshaw , S. & Fischer , A. J. Midkine is neuroprotective and influences glial reactivity and the formation Müller glia-derived progenitor cells in chick and mouse retinas . Glia 69 , 1515 – 1539 ( 2021 ). OpenUrl CrossRef PubMed 88. ↵ Murphy , M. P. & LeVine , H. Alzheimer’s Disease and the β-Amyloid Peptide . J. Alzheimers Dis. JAD 19 , 311 ( 2010 ). OpenUrl PubMed 89. ↵ Metaxas , A. & Kempf , S. J. Neurofibrillary tangles in Alzheimer’s disease: elucidation of the molecular mechanism by immunohistochemistry and tau protein phospho-proteomics . Neural Regen. Res . 11 , 1579 – 1581 ( 2016 ). OpenUrl PubMed 90. ↵ Keren-Shaul , H. , Spinrad , A. , Weiner , A. , Matcovitch-Natan , O. , Dvir-Szternfeld , R. , Ulland , T. K. , David , E. , Baruch , K. , Lara-Astaiso , D. , Toth , B. , Itzkovitz , S. , Colonna , M. , Schwartz , M. & Amit , I. A Unique Microglia Type Associated with Restricting Development of Alzheimer’s Disease . Cell 169 , 1276 – 1290.e17 ( 2017 ). OpenUrl CrossRef PubMed 91. ↵ Leng , F. & Edison , P. Neuroinflammation and microglial activation in Alzheimer disease: where do we go from here? Nat. Rev. Neurol . 17 , 157 – 172 ( 2021 ). OpenUrl CrossRef PubMed 92. ↵ Hua , H. , Münter , L. , Harmeier , A. , Georgiev , O. , Multhaup , G. & Schaffner , W. Toxicity of Alzheimer’s disease-associated Aβ peptide is ameliorated in a Drosophila model by tight control of zinc and copper availability . Biol. Chem . 392 , 919 – 926 ( 2011 ). OpenUrl PubMed 93. ↵ Zhang , Z. , Deng , X. , Liu , Y. , Liu , Y. , Sun , L. & Chen , F. PKM2, function and expression and regulation . Cell Biosci . 9 , 52 ( 2019 ). OpenUrl CrossRef PubMed 94. ↵ Cao , X. , Sun , H. , Feng , R. , Mazumder , R. , Najar , C. F. B. A. , Li , Y. I. , Jager, P. L. de , Bennett , D. , Consortium , T. A. D. F. G. , Dey , K. K. & Wang , G. Integrative multi-omics QTL colocalization maps regulatory architecture in aging human brain . 2025.04.17.25326042 Preprint at doi: 10.1101/2025.04.17.25326042 ( 2025 ) OpenUrl Abstract / FREE Full Text 95. ↵ Kazmi , S. Z. H. & Wierenga , C. J. Deep cortical layers are more vulnerable in Alzheimer’s disease . Neuron 113 , 2218 – 2220 ( 2025 ). OpenUrl PubMed 96. ↵ Braak , H. & Braak , E. Neuropathological stageing of Alzheimer-related changes . Acta Neuropathol. (Berl .) 82 , 239 – 259 ( 1991 ). OpenUrl CrossRef PubMed Web of Science 97. ↵ Forner , S. , Kawauchi , S. , Balderrama-Gutierrez , G. , Kramár , E. A. , Matheos , D. P. , Phan , J. , Javonillo , D. I. , Tran , K. M. , Hingco , E. , da Cunha , C. , Rezaie , N. , Alcantara , J. A. , Baglietto-Vargas , D. , Jansen , C. , Neumann , J. , Wood , M. A. , MacGregor , G. R. , Mortazavi , A. , Tenner , A. J. , LaFerla , F. M. & Green , K. N. Systematic phenotyping and characterization of the 5xFAD mouse model of Alzheimer’s disease . Sci. Data 8 , 270 ( 2021 ). OpenUrl PubMed 98. ↵ Lish , A. M. , Grogan , E. F. L. , Benoit , C. R. , Pearse , R. V. , Heuer , S. E. , Luquez , T. , Orme , G. A. , Galle , P. C. , Milinkeviciute , G. , Green , K. N. , Alexander , K. D. , Fancher , S. B. , Stern , A. M. , Fujita , M. , Bennett , D. A. , Seyfried , N. T. , Jager , P. L. D. , Menon , V. & Young-Pearse , T. L. CLU alleviates Alzheimer’s disease-relevant processes by modulating astrocyte reactivity and microglia-dependent synaptic density . Neuron 113 , 1925 – 1946.e11 ( 2025 ). OpenUrl PubMed 99. ↵ Foster , E. M. , Dangla-Valls , A. , Lovestone , S. , Ribe , E. M. & Buckley , N. J. Clusterin in Alzheimer’s Disease: Mechanisms, Genetics, and Lessons From Other Pathologies . Front. Neurosci . 13 , 164 ( 2019 ). OpenUrl CrossRef PubMed 100. ↵ Fu , C.-Y. , Chen , H.-Y. , Lin , C.-Y. , Chen , S.-J. Sheu , J.-C. & Tsai , H.-J. Extracellular Pgk1 interacts neural membrane protein enolase-2 to improve the neurite outgrowth of motor neurons . Commun. Biol . 6 , 849 ( 2023 ). OpenUrl PubMed 101. ↵ Zhang , K. , Sun , L. & Kang , Y. Regulation of phosphoglycerate kinase 1 and its critical role in cancer . Cell Commun. Signal. CCS 21 , 240 ( 2023 ). OpenUrl PubMed 102. ↵ Wang , J. , Li , J. , Kramer , S. T. , Su , L. , Chang , Y. , Xu , C. , Eadon , M. T. , Kiryluk , K. , Ma , Q. & Xu , D. Dimension-agnostic and granularity-based spatially variable gene identification using BSP . Nat. Commun . 14 , 7367 ( 2023 ). OpenUrl CrossRef PubMed 103. ↵ Liu , Q. Hsu , C.-Y. & Shyr , Y. Scalable and model-free detection of spatial patterns and colocalization . Genome Res . 32 , 1736 – 1745 ( 2022 ). OpenUrl Abstract / FREE Full Text 104. ↵ Chen , Z. , Soifer , I. , Hilton , H. , Keren , L. & Jojic , V. Modeling Multiplexed Images with Spatial-LDA Reveals Novel Tissue Microenvironments . J. Comput. Biol . 27 , 1204 – 1218 ( 2020 ). OpenUrl CrossRef PubMed 105. ↵ Miller , B. F. , Huang , F. , Atta , L. , Sahoo , A. & Fan , J. Reference-free cell type deconvolution of multi-cellular pixel-resolution spatially resolved transcriptomics data . Nat. Commun . 13 , 2339 ( 2022 ). OpenUrl CrossRef PubMed 106. ↵ Cang , Z. , Zhao , Y. , Almet , A. A. , Stabell , A. , Ramos , R. , Plikus , M. V. , Atwood , S. X. & Nie , Q. Screening cell–cell communication in spatial transcriptomics via collective optimal transport . Nat. Methods 20 , 218 – 228 ( 2023 ). OpenUrl CrossRef PubMed 107. ↵ Browaeys , R. , Saelens , W. & Saeys , Y. NicheNet: modeling intercellular communication by linking ligands to target genes . Nat. Methods 17 , 159 – 162 ( 2020 ). OpenUrl CrossRef PubMed 108. ↵ Fang , H. , ULTRA-DD Consortium , De Wolf , H. , Knezevic , B. , Burnham , K. L. , Osgood , J. , Sanniti , A. , Lledó Lara , A. , Kasela , S. , De Cesco , S. , Wegner , J. K. , Handunnetthi , L. , McCann , F. E. , Chen , L. , Sekine , T. , Brennan , P. E. , Marsden , B. D. , Damerell , D. , O’Callaghan , C. A. , Bountra , C. , Bowness , P. , Sundström , Y. , Milani , L. , Berg , L. , Göhlmann , H. W. , Peeters , P. J. , Fairfax , B. P. , Sundström , M. & Knight , J. C. A genetics-led approach defines the drug target landscape of 30 immune-related traits . Nat. Genet . 51 , 1082 – 1091 ( 2019 ). OpenUrl CrossRef PubMed 109. ↵ Mancuso , N. , Freund , M. K. , Johnson , R. , Shi , H. , Kichaev , G. , Gusev , A. & Pasaniuc , B. Probabilistic fine-mapping of transcriptome-wide association studies . Nat. Genet . 51 , 675 – 682 ( 2019 ). OpenUrl CrossRef PubMed 110. ↵ DeTomaso , D. & Yosef , N. Hotspot identifies informative gene modules across modalities of single-cell genomics . Cell Syst . 12 , 446 – 456.e9 ( 2021 ). OpenUrl CrossRef PubMed 111. ↵ Efremova , M. , Vento-Tormo , M. , Teichmann , S. A. & Vento-Tormo , R. CellPhoneDB: inferring cell– cell communication from combined expression of multi-subunit ligand–receptor complexes . Nat. Protoc . 15 , 1484 – 1506 ( 2020 ). OpenUrl CrossRef PubMed 112. ↵ Xuanyuan , Q. , Wu , H. , Sundaramoorthi , H. , Isnard , P. , Chen , C. , Rahmani , W. & Humphreys , B. D. Multimodal spatial transcriptomic characterization of mouse kidney injury and repair . Nat. Commun . 16 , 7567 ( 2025 ). OpenUrl PubMed 113. ↵ Janesick , A. , Shelansky , R. , Gottscho , A. D. , Wagner , F. , Williams , S. R. , Rouault , M. , Beliakoff , G. , Morrison , C. A. , Oliveira , M. F. , Sicherman , J. T. , Kohlway , A. , Abousoud , J. , Drennon , T. Y. , Mohabbat , S. H. & Taylor , S. E. B. High resolution mapping of the tumor microenvironment using integrated single-cell, spatial and in situ analysis . Nat. Commun . 14 , 8353 ( 2023 ). OpenUrl CrossRef PubMed 114. ↵ Zhao , P. , Zhu , J. , Ma , Y. & Zhou , X. Modeling zero inflation is not necessary for spatial transcriptomics . Genome Biol . 23 , 118 ( 2022 ). OpenUrl CrossRef PubMed 115. ↵ Shang , L. , Wu , P. & Zhou , X. Statistical identification of cell type-specific spatially variable genes in spatial transcriptomics . Nat. Commun . 16 , 1059 ( 2025 ). OpenUrl CrossRef PubMed 116. ↵ Zhang , D. & Lin , X. Hypothesis testing in semiparametric additive mixed models . Biostatistics 4 , 57 – 74 ( 2003 ). OpenUrl CrossRef PubMed Web of Science 117. ↵ 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 ( 2020 ). OpenUrl CrossRef PubMed 118. ↵ Bhuva , D. D. , Tan , C. W. , Salim , A. , Marceaux , C. , Pickering , M. A. , Chen , J. , Kharbanda , M. , Jin , X. , Liu , N. , Feher , K. , Putri , G. , Tilley , W. D. , Hickey , T. E. , Asselin-Labat , M.-L. , Phipson , B. & Davis , M. J. Library size confounds biology in spatial transcriptomics data . Genome Biol . 25 , 99 ( 2024 ). OpenUrl CrossRef PubMed 119. ↵ Salim , A. , Bhuva , D. D. , Chen , C. , Tan , C. W. , Yang , P. , Davis , M. J. & Yang , J. Y. H. SpaNorm: spatially-aware normalization for spatial transcriptomics data . Genome Biol . 26 , 109 ( 2025 ). OpenUrl PubMed 120. ↵ McCarthy , D. J. , Campbell , K. R. , Lun , A. T. L. & Wills , Q. F. Scater: pre-processing, quality control, normalization and visualization of single-cell RNA-seq data in R . Bioinformatics 33 , 1179 – 1186 ( 2017 ). OpenUrl CrossRef PubMed 121. ↵ Palla , G. , Spitzer , H. , Klein , M. , Fischer , D. , Schaar , A. C. , Kuemmerle , L. B. , Rybakov , S. , Ibarra , I. L. , Holmberg , O. , Virshup , I. , Lotfollahi , M. , Richter , S. & Theis , F. J. Squidpy: a scalable framework for spatial omics analysis . Nat. Methods 19 , 171 – 178 ( 2022 ). OpenUrl CrossRef PubMed 122. ↵ Svensson , V. , Teichmann , S. A. & Stegle , O. SpatialDE: identification of spatially variable genes . Nat. Methods 15 , 343 – 346 ( 2018 ). OpenUrl CrossRef PubMed 123. ↵ Kats , I. , Vento-Tormo , R. & Stegle , O. SpatialDE2: Fast and localized variance component analysis of spatial transcriptomics . Preprint at doi: 10.1101/2021.10.27.466045 ( 2021 ) OpenUrl Abstract / FREE Full Text 124. ↵ Weber , L. M. , Saha , A. , Datta , A. , Hansen , K. D. & Hicks , S. C. nnSVG for the scalable identification of spatially variable genes using nearest-neighbor Gaussian processes . Nat. Commun . 14 , 4059 ( 2023 ). OpenUrl CrossRef PubMed 125. ↵ Sun , S. , Zhu , J. & Zhou , X. Statistical Analysis of Spatial Expression Pattern for Spatially Resolved Transcriptomic Studies . Preprint at doi: 10.1101/810903 ( 2019 ) OpenUrl Abstract / FREE Full Text 126. ↵ Zhu , J. , Sun , S. & Zhou , X. SPARK-X: non-parametric modeling enables scalable and robust detection of spatial expression patterns for large spatial transcriptomic studies . Genome Biol . 22 , 184 ( 2021 ). OpenUrl CrossRef PubMed 127. ↵ Su , H. , Wu , Y. , Chen , B. & Cui , Y. STANCE: a unified statistical model to detect cell-type-specific spatially variable genes in spatial transcriptomics . Preprint at doi: 10.1101/2024.09.22.614385 ( 2024 ) OpenUrl Abstract / FREE Full Text 128. ↵ Song , L. , Chen , W. , Hou , J. , Guo , M. & Yang , J. Spatially resolved mapping of cells associated with human complex traits . Preprint at doi: 10.1101/2024.10.31.24316538 ( 2024 ) OpenUrl Abstract / FREE Full Text 129. ↵ de Leeuw , C. A. , Mooij , J. M. , Heskes , T. & Posthuma , D. MAGMA: Generalized Gene-Set Analysis of GWAS Data . PLoS Comput. Biol . 11 , e1004219 ( 2015 ). OpenUrl CrossRef PubMed 130. ↵ Weeks , E. M. , Ulirsch , J. C. , Cheng , N. Y. , Trippe , B. L. , Fine , R. S. , Miao , J. , Patwardhan , T. A. , Kanai , M. , Nasser , J. , Fulco , C. P. , Tashman , K. C. , Aguet , F. , Li , T. , Ordovas-Montanes , J. , Smillie , C. S. , Biton , M. , Shalek , A. K. , Ananthakrishnan , A. N. , Xavier , R. J. , Regev , A. , Gupta , R. M. , Lage , K. , Ardlie , K. G. , Hirschhorn , J. N. , Lander , E. S. , Engreitz , J. M. & Finucane , H. K. Leveraging polygenic enrichments of gene features to predict genes underlying complex traits and diseases . ( 2024 ). 131. ↵ L. Lun , A.T., Bach , K. & Marioni , J. C. Pooling across cells to normalize single-cell RNA sequencing data with many zero counts . Genome Biol . 17 , 75 ( 2016 ). OpenUrl CrossRef PubMed 132. ↵ Kim , S. S. , Dai , C. , Hormozdiari , F. , Geijn, B. van de , Gazal , S. , Park , Y. , O’Connor , L. , Amariuta , T. , Loh , P.-R. , Finucane , H. , Raychaudhuri , S. & Price , A. L. Genes with High Network Connectivity Are Enriched for Disease Heritability . Am. J. Hum. Genet . 104 , 896 – 913 ( 2019 ). OpenUrl CrossRef PubMed 133. ↵ Cao , X. , Sun , H. , Feng , R. , Mazumder , R. , Najar , C. F. B. A. , Li , Y. I. , Jager , P.L.de , Bennett , D. , Consortium , T. A. D. F. G. , Dey , K. K. & Wang , G. Integrative multi-omics QTL colocalization maps regulatory architecture in aging human brain . 2025.04.17.25326042 Preprint at doi: 10.1101/2025.04.17.25326042 ( 2025 ) OpenUrl Abstract / FREE Full Text 134. ↵ Najar , C.F.B.A. , Feng , R. , Dai , C. , Fair , B. , Hauck , Q. , Li , J. , Cao , X. , Dey , K.K. , De Jager , P. , Bennett , D. , Liu , X. , Wang , G. & Li , Y.I. Genetic and functional analysis of unproductive splicing using LeafCutter2 . bioRxiv 2025.04.06.646893 ( 2025 ). doi: 10.1101/2025.04.06.646893 OpenUrl Abstract / FREE Full Text 135. ↵ Bellenguez , C. , Küçükali , F. , Jansen , I. E. , Kleineidam , L. , Moreno-Grau , S. , Amin , N. , Naj , A. C. , Campos-Martin , R. , Grenier-Boley , B. , Andrade , V. , Holmans , P. A. , Boland , A. , Damotte , V. , van der Lee , S. J. , Costa , M. R. , Kuulasmaa , T. , Yang , Q. , de Rojas , I. , Bis , J. C. , Yaqub , A. , Prokic , I. , Chapuis , J. , Ahmad , S. , Giedraitis , V. , Aarsland , D. , Garcia-Gonzalez , P. , Abdelnour , C. , Alarcón-Martín , E. , Alcolea , D. , Alegret , M. , Alvarez , I. , Álvarez , V. , Armstrong , N. J. , Tsolaki , A. , Antúnez , C. , Appollonio , I. , Arcaro , M. , Archetti , S. , Pastor , A. A. , Arosio , B. , Athanasiu , L. , Bailly , H. , Banaj , N. , Baquero , M. , Barral , S. , Beiser , A. , Pastor , A. B. , Below , J. E. , Benchek , P. , Benussi , L. , Berr , C. , Besse , C. , Bessi , V. , Binetti , G. , Bizarro , A. , Blesa , R. , Boada , M. , Boerwinkle , E. , Borroni , B. , Boschi , S. , Bossù , P. , Bråthen , G. , Bressler , J. , Bresner , C. , Brodaty , H. , Brookes , K. J. , Brusco , L. I. , Buiza-Rueda , D. , Bûrger , K. , Burholt , V. , Bush , W. S. , Calero , M. , Cantwell , L. B. , Chene , G. , Chung , J. , Cuccaro , M. L. , Carracedo , Á. , Cecchetti , R. , Cervera-Carles , L. , Charbonnier , C. , Chen , H.-H. , Chillotti , C. , Ciccone , S. , Claassen , J. A. H. R. , Clark , C. , Conti , E. , Corma-Gómez , A. , Costantini , E. , Custodero , C. , Daian , D. , Dalmasso , M. C. , Daniele , A. , Dardiotis , E. , Dartigues , J.-F. , de Deyn , P. P. , de Paiva Lopes , K. , de Witte , L. D. , Debette , S. , Deckert , J. , del Ser , T. , Denning , N. , DeStefano , A. , Dichgans , M. , Diehl-Schmid , J. , Diez-Fairen , M. , Rossi , P. D. , Djurovic , S. , Duron , E. , Düzel , E. , Dufouil , C. , Eiriksdottir , G. , Engelborghs , S. , Escott-Price , V. , Espinosa , A. , Ewers , M. , Faber , K. M. , Fabrizio , T. , Nielsen , S. F. , Fardo , D. W. , Farotti , L. , Fenoglio , C. , Fernández-Fuertes , M. , Ferrari , R. , Ferreira , C. B. , Ferri , E. , Fin , B. , Fischer , P. , Fladby , T. , Fließbach , K. , Fongang , B. , Fornage , M. , Fortea , J. , Foroud , T. M. , Fostinelli , S. , Fox , N. C. , Franco-Macías , E. , Bullido , M. J. , Frank-García , A. , Froelich , L. , Fulton-Howard , B. , Galimberti , D. , García-Alberca , J. M. , García-González , P. , Garcia-Madrona , S. , Garcia-Ribas , G. , Ghidoni , R. , Giegling , I. , Giorgio , G. , Goate , A. M. , Goldhardt , O. , Gomez-Fonseca , D. , González-Pérez , A. , Graff , C. , Grande , G. , Green , E. , Grimmer , T. , Grünblatt , E. , Grunin , M. , Gudnason , V. , Guetta-Baranes , T. , Haapasalo , A. , Hadjigeorgiou , G. , Haines , J. L. , Hamilton-Nelson , K. L. , Hampel , H. , Hanon , O. , Hardy , J. , Hartmann , A. M. , Hausner , L. , Harwood , J. , Heilmann-Heimbach , S. , Helisalmi , S. , Heneka , M. T. , Hernández , I. , Herrmann , M. J. , Hoffmann , P. , Holmes , C. , Holstege , H. , Vilas , R. H. , Hulsman , M. , Humphrey , J. , Biessels , G. J. , Jian , X. , Johansson , C. , Jun , G. R. , Kastumata , Y. , Kauwe , J. , Kehoe , P. G. , Kilander , L. , Ståhlbom , A. K. , Kivipelto , M. , Koivisto , A. , Kornhuber , J. , Kosmidis , M. H. , Kukull , W. A. , Kuksa , P. P. , Kunkle , B. W. , Kuzma , A. B. , Lage , C. , Laukka , E. J. , Launer , L. , Lauria , A. , Lee , C.-Y. , Lehtisalo , J. , Lerch , O. , Lleó , A. , Longstreth , W. , Lopez , O. , de Munain , A. L. , Love , S. , Löwemark , M. , Luckcuck , L. , Lunetta , K. L. , Ma , Y. , Macías , J. , MacLeod , C. A. , Maier , W. , Mangialasche , F. , Spallazzi , M. , Marquié , M. , Marshall , R. , Martin , E. R. , Montes , A. M. , Rodríguez , C. M. , Masullo , C. , Mayeux , R. , Mead , S. , Mecocci , P. , Medina , M. , Meggy , A. , Mehrabian , S. , Mendoza , S. , Menéndez-González , M. , Mir , P. , Moebus , S. , Mol , M. , Molina-Porcel , L. , Montrreal , L. , Morelli , L. , Moreno , F. , Morgan , K. , Mosley , T. , Nöthen , M. M. , Muchnik , C. , Mukherjee , S. , Nacmias , B. , Ngandu , T. , Nicolas , G. , Nordestgaard , B. G. , Olaso , R. , Orellana , A. , Orsini , M. , Ortega , G. , Padovani , A. , Paolo , C. , Papenberg , G. , Parnetti , L. , Pasquier , F. , Pastor , P. , Peloso , G. , Pérez-Cordón , A. , Pérez-Tur , J. , Pericard , P. , Peters , O. , Pijnenburg , Y. A. L. , Pineda , J. A. , Piñol-Ripoll , G. , Pisanu , C. , Polak , T. , Popp , J. , Posthuma , D. , Priller , J. , Puerta , R. , Quenez , O. , Quintela , I. , Thomassen , J. Q. , Rábano , A. , Rainero , I. , Rajabli , F. , Ramakers , I. , Real , L. M. , Reinders , M. J. T. , Reitz , C. , Reyes-Dumeyer , D. , Ridge , P. , Riedel-Heller , S. , Riederer , P. , Roberto , N. , Rodriguez-Rodriguez , E. , Rongve , A. , Allende , I. R. , Rosende-Roca , M. , Royo , J. L. , Rubino , E. , Rujescu , D. , Sáez , M. E. , Sakka , P. , Saltvedt , I. , Sanabria , Á. , Sánchez-Arjona , M. B. , Sanchez-Garcia , F. , Juan , P. S. , Sánchez-Valle , R. , Sando , S. B. , Sarnowski , C. , Satizabal , C. L. , Scamosci , M. , Scarmeas , N. , Scarpini , E. , Scheltens , P. , Scherbaum , N. , Scherer , M. , Schmid , M. , Schneider , A. , Schott , J. M. , Selbæk , G. , Seripa , D. , Serrano , M. , Sha , J. , Shadrin , A. A. , Skrobot , O. , Slifer , S. , Snijders , G. J. L. , Soininen , H. , Solfrizzi , V. , Solomon , A. , Song , Y. , Sorbi , S. , Sotolongo-Grau , O. , Spalletta , G. , Spottke , A. , Squassina , A. , Stordal , E. , Tartan , J. P. , Tárraga , L. , Tesí , N. , Thalamuthu , A. , Thomas , T. , Tosto , G. , Traykov , L. , Tremolizzo , L. , Tybjærg-Hansen , A. , Uitterlinden , A. , Ullgren , A. , Ulstein , I. , Valero , S. , Valladares , O. , Broeckhoven , C. V. , Vance , J. , Vardarajan , B. N. , van der Lugt , A. , Dongen , J. V. , van Rooij , J. , van Swieten , J. , Vandenberghe , R. , Verhey , F. , Vidal , J.-S. , Vogelgsang , J. , Vyhnalek , M. , Wagner , M. , Wallon , D. , Wang , L.-S. , Wang , R. , Weinhold , L. , Wiltfang , J. , Windle , G. , Woods , B. , Yannakoulia , M. , Zare , H. , Zhao , Y. , Zhang , X. , Zhu , C. , Zulaica , M. , Farrer , L. A. , Psaty , B. M. , Ghanbari , M. , Raj , T. , Sachdev , P. , Mather , K. , Jessen , F. , Ikram , M. A. , de Mendonça , A. , Hort , J. , Tsolaki , M. , Pericak-Vance , M. A. , Amouyel , P. , Williams , J. , Frikke-Schmidt , R. , Clarimon , J. , Deleuze , J.-F. , Rossi , G. , Seshadri , S. , Andreassen , O. A. , Ingelsson , M. , Hiltunen , M. , Sleegers , K. , Schellenberg , G. D. , van Duijn , C. M. , Sims , R. , van der Flier , W. M. , Ruiz , A. , Ramirez , A. & Lambert , J.-C. New insights into the genetic etiology of Alzheimer’s disease and related dementias . Nat. Genet . 54 , 412 – 436 ( 2022 ). OpenUrl CrossRef PubMed 136. ↵ Jin , X. , Simmons , S. K. , Guo , A. , Shetty , A. S. , Ko , M. , Nguyen , L. , Jokhi , V. , Robinson , E. , Oyler , P. , Curry , N. , Deangeli , G. , Lodato , S. , Levin , J. Z. , Regev , A. , Zhang , F. & Arlotta , P. In vivo Perturb-Seq reveals neuronal and glial abnormalities associated with autism risk genes . Science 370 , eaaz6063 ( 2020 ). OpenUrl Abstract / FREE Full Text View the discussion thread. Back to top Previous Next Posted September 25, 2025. Download PDF Supplementary Material Email Thank you for your interest in spreading the word about bioRxiv. NOTE: Your email address is requested solely to identify you as the sender of this article. Your Email * Your Name * Send To * Enter multiple addresses on separate lines or separate them with commas. You are going to email the following Mapping disease critical spatially variable gene programs by integrating spatial transcriptomics with human genetics Message Subject (Your Name) has forwarded a page to you from bioRxiv Message Body (Your Name) thought you would like to see this page from the bioRxiv website. Your Personal Message CAPTCHA This question is for testing whether or not you are a human visitor and to prevent automated spam submissions. Share Mapping disease critical spatially variable gene programs by integrating spatial transcriptomics with human genetics Hanbyul Lee , Haochen Sun , Xuewei Cao , Berke Karaahmet , Zhijian Li , Hans Ulrich-Klein , Mariko Taga , Gao Wang , Philip L. De Jager , David A. Bennett , Luca Pinello , Xin Jin , Rahul Mazumder , Kushal K. Dey bioRxiv 2025.09.24.678397; doi: https://doi.org/10.1101/2025.09.24.678397 Share This Article: Copy Citation Tools Mapping disease critical spatially variable gene programs by integrating spatial transcriptomics with human genetics Hanbyul Lee , Haochen Sun , Xuewei Cao , Berke Karaahmet , Zhijian Li , Hans Ulrich-Klein , Mariko Taga , Gao Wang , Philip L. De Jager , David A. Bennett , Luca Pinello , Xin Jin , Rahul Mazumder , Kushal K. Dey bioRxiv 2025.09.24.678397; doi: https://doi.org/10.1101/2025.09.24.678397 Citation Manager Formats BibTeX Bookends EasyBib EndNote (tagged) EndNote 8 (xml) Medlars Mendeley Papers RefWorks Tagged Ref Manager RIS Zotero Tweet Widget Facebook Like Google Plus One Subject Area Genomics Subject Areas All Articles Animal Behavior and Cognition (7637) Biochemistry (17705) Bioengineering (13899) Bioinformatics (41968) Biophysics (21460) Cancer Biology (18603) Cell Biology (25526) Clinical Trials (138) Developmental Biology (13385) Ecology (19909) Epidemiology (2067) Evolutionary Biology (24326) Genetics (15614) Genomics (22513) Immunology (17741) Microbiology (40423) Molecular Biology (17193) Neuroscience (88645) Paleontology (667) Pathology (2835) Pharmacology and Toxicology (4825) Physiology (7647) Plant Biology (15160) Scientific Communication and Education (2046) Synthetic Biology (4302) Systems Biology (9825) Zoology (2271)