Coarsening Visium HD resolution from 8 to 64 μm can flip cell-type colocalization from negative to positive ( r = - 0.12 → + 0.80 ) , yet many widely used compositional deconvolution workflows require coarsening or subsampling at million-bin scale. Here we introduce FlashDeconv, which combines leverage-score importance sampling with sparse spatial regularization to achieve competitive benchmark accuracy while processing 1.6 million bins in 153 seconds on commodity hardware. Systematic multi-resolution analysis of Visium HD mouse intestine reveals a tissue-specific resolution horizon (8-16 μ m)-the scale at which this sign inversion occurs-validated by Xenium ground truth. Below this horizon, FlashDeconv provides, to our knowledge, the first sequencing-based quantification of Tuft cell chemosensory niches (15.3-fold stem cell enrichment). In a 1.6-million-bin human colorectal cancer cohort, FlashDeconv uncovers neutrophil inflammatory microdomains co-localized with immunoregulatory dendritic cells (mRegDC) at the tumor-stroma interface-spatial niches largely missed by discrete-label summaries, with RCTD doublet mode labeling only 2.3% of hotspot bins as neutrophil singlets.

The rapid expansion of single-cell RNA sequencing (scRNA-seq) has made accurate cell type annotation a critical bottleneck for biological discovery. Existing computational methods are often limited by reference data dependency, while emerging single Large Language Model (LLM) approaches are susceptible to model-specific biases and provide insufficient uncertainty quantification. To address these limitations, we introduce mLLMCelltype, a framework that harnesses collective intelligence-the emergent problem-solving capacity arising when multiple independent agents interact through structured deliberation to produce solutions exceeding individual capabilities-of multiple LLMs through an iterative deliberation process. Across 49 diverse datasets, our framework achieves a mean accuracy of 77.2%, a 15.7-percentage-point improvement over the best-performing single-LLM baseline (61.5%). The consensus mechanism demonstrates high robustness to noisy input and generalizes to datasets released after the LLMs' training. By providing transparent reasoning chains and robust consensus-based confidence metrics, mLLMCelltype minimizes manual annotation effort and enables reliable interpretation of complex cellular landscapes. The framework is available as an open-source package and an accessible web server.

Spatial transcriptomics has transformed our ability to study tissue architecture at molecular resolution, yet analyzing these data demands navigating dozens of computational methods across incompatible Python and R ecosystems-forcing researchers to devote more effort to making tools function than to pursuing biological questions. We present ChatSpatial, a platform in which the LLM selects from pre-validated tool schemas rather than generating free-form code, with domain expertise embedded in schema descriptions for context-aware parameter inference. Built on the Model Context Protocol (MCP), ChatSpatial unifies 60+ methods across 15 analytical categories into a single conversational workflow spanning Python and R ecosystems. Replication of two published studies-recovering subclonal heterogeneity in ovarian cancer and tumor microenvironment organization in oral squamous cell carcinoma-and validation across seven LLM platforms demonstrate that schema-enforced orchestration yields near-deterministic reproducibility at the workflow level for multi-step spatial analyses. Beyond replication, exploratory cross-method analyses illustrate practical triangulation across independent analytical frameworks.

In genome-wide epigenetic studies, determining how exposures (e.g., Single Nucleotide Polymorphisms) affect outcomes (e.g., gene expression) through intermediate variables, such as DNA methylation, is a key challenge. Mediation analysis provides a framework to identify these causal pathways; however, testing for mediation effects involves a complex composite null hypothesis. Existing methods, such as Sobel's test or the Max-P test, are often underpowered in this context because they rely on null distributions determined under only a subset of the null space and are not optimized for the multiple testing burden inherent in high-dimensional data. To address these limitations, we introduce MLFDR (Mediation Analysis using Local False Discovery Rates), a novel method for high-dimensional mediation analysis. MLFDR leverages local false discovery rates, calculated from the coefficients of structural equation models, to construct an optimal rejection region. We demonstrate theoretically and through simulation that MLFDR asymptotically controls the false discovery rate and achieves superior statistical power compared to recent high-dimensional mediation methods. In real data applications, MLFDR identified 20%-50% more significant mediators than existing methods, demonstrating its ability to uncover biological signals missed by conventional approaches.

Methods that map genetic risk to cells do not directly test whether spatially organized ligand-receptor (LR) gene annotations carry conditional heritability association. Here we introduce EdgeMap, which scores each gene by the spatial activity of its LR contexts in spatial transcriptomics data, maps these scores alongside cell-intrinsic annotations to SNP-level LD scores, and tests both jointly against GWAS summary statistics. Across 17 traits and five human tissues, edge z -scores are systematically higher in biologically matched tissues: all 15 traits with both matched and unmatched tissue classes show this ordering (median Δ z = 1.51 ; sign P = 3.1 × 10 -5 ), and 11 of 13 nominal-positive associations concentrate in the biologically matched set ( P = 1 × 10 -4 ). The pattern survives broad LR-gene-class and cell-type composition controls. Donor-level cross-section summaries, independent GWAS replication for LDL and CAD, and cell-segmented Visium HD liver data provide convergent support. A secondary analysis prioritizes constituent genes within active LR contexts; sixteen of 31 prioritized genes are absent from standard gene-level methods. These results establish spatially weighted LR-gene annotations as a complementary layer for interpreting complex-trait genetic architecture.

Identifying spatially variable genes (SVGs) is the first analytical step in spatial transcriptomics, determining which genes and pathways are prioritized for downstream validation. Yet the restricted spatial models of current detection methods create systematic blind spots that can exclude biologically coherent programs from discovery. Here we present FlashS, which reformulates kernel-based spatial testing in the frequency domain to detect arbitrary multi-scale expression patterns while scaling to millions of cells. In human cardiac tissue, this broader detection capacity recovers a coherent PGC-1α-regulated mitochondrial biogenesis program-40 of 49 pathway genes spatially associated with ventricular cardiomyocytes-that PreTSA, a leading parametric alternative, largely misses (1 of 49 genes), a finding replicated in an independent cohort. Across 50 benchmark datasets spanning 9 platforms, FlashS achieves state-of-the-art ranking accuracy (mean Kendall τ = 0.935 ) and completes on the Allen Brain MERFISH atlas (3.94 million cells) in 12.6 minutes with 21.5 GB memory.

Many multiple-testing methods compare the two sides of a null distribution to control the false discovery rate (FDR). Small p-values or large positive scores are treated as possible discoveries, while large p-values or large negative scores are used to estimate how many of those discoveries are false. The mirror and knockoff+ thresholds are built on this idea. For valid knockoff statistics, the comparison is justified by a strong property: conditional on their magnitudes and the nonnull scores, the null signs are independent fair coin flips. This paper asks what can go wrong when the same threshold is used without that property. We give three answers. First, we construct exactly uniform p-values satisfying positive regression dependence on a subset (PRDS), with a joint density that is positive throughout the unit cube. At a nominal level of 10 These results do not contradict knockoff theory. They show that the shared counting rule is not, by itself, an FDR guarantee: validity depends on the joint behavior of the null signs, not only on marginal symmetry, Gaussianity, or PRDS.

Generative, temporal network models play an important role in analyzing the dependence structure and evolution patterns of complex networks. Due to the complicated nature of real network data, it is often naive to assume that the underlying data-generative mechanism itself is invariant with time. Such observation leads to the study of changepoints or sudden shifts in the distributional structure of the evolving network. In this paper, we propose a likelihood-based methodology to detect changepoints in undirected, affine preferential attachment networks, and establish a hypothesis testing framework to detect a single changepoint, together with a consistent estimator for the changepoint. Such results require establishing consistency and asymptotic normality of the MLE under the changepoint regime, which suffers from long range dependence. The methodology is then extended to the multiple changepoint setting via both a sliding window method and a more computationally efficient score statistic. We also compare the proposed methodology with previously developed non-parametric estimators of the changepoint via simulation, and the methods developed herein are applied to modeling the popularity of a topic in a Twitter network over time.

The unbiased sample versions of squared distance covariance and the Hilbert-Schmidt independence criterion (HSIC) are fourth-order U-statistics, yet U-centering evaluates them from pairwise arrays in O(n^2) operations. We show that U-centering is exactly the least-squares residual obtained after fitting additive endpoint effects to a symmetric hollow array. This interpretation explains the zero row sums and the denominator n(n-3) through the residual degrees of freedom. The same pairwise residualization also gives useful regression identities. After endpoint effects are removed from both arrays, the U-centered dependence t-statistic is the ordinary slope t-statistic obtained by regressing one adjusted array on the other. In the two-sample problem, pooling the observations and using the between-group pair indicator as the predictor shows that the generalized-energy statistic is twice the fitted slope. The common-endpoint and fully interacted regressions give the same slope but use different residual standard errors. For n≥2r, we extend the construction to arrays indexed by r-subsets. Higher-order U-centering removes all effects involving fewer than r sample labels, leaves zero (r-1)-way margins, and projects onto a residual space of dimension nr- nr-1. For two symmetric kernels with r arguments, the normalized inner product of the centered arrays is unbiased for the cross-moment of their rth Hoeffding components. A direct estimator can involve products spanning as many as 2r observations, but subset-margin inversion or higher-order U-centering evaluates the same quantity in O(n^r) operations for fixed r. When both arrays are formed from the same kernel and sample, this becomes a nonnegative unbiased estimator of the variance of the highest-order Hoeffding component.

Change point analysis is concerned with detecting and locating structure breaks in the underlying model of a sequence of observations ordered by time, space or other variables. A widely adopted approach for change point analysis is to minimize an objective function with a penalty term on the number of change points. This framework includes several well-established procedures, such as the penalized log-likelihood using the (modified) Bayesian information criterion (BIC) or the minimum description length (MDL). The resulting optimization problem can be solved in polynomial time by dynamic programming or its improved version, such as the Pruned Exact Linear Time (PELT) algorithm (Killick, Fearnhead, and Eckley 2012). However, existing computational methods often suffer from two primary limitations: (1) methods based on direct implementation of dynamic programming or PELT are often time-consuming for long data sequences due to repeated computation of the cost value over different segments of the data sequence; (2) state-of-the-art R packages do not provide enough flexibility for users to handle different change point settings and models. In this work, we present the fastcpd package, aiming to provide an efficient and versatile framework for change point detection in several commonly encountered settings. The core of our algorithm is built upon PELT and the sequential gradient descent method recently proposed by Zhang and Dawn (2023). We illustrate the usage of the fastcpd package through several examples, including mean/variance changes in a (multivariate) Gaussian sequence, parameter changes in regression models, structural breaks in ARMA/GARCH/VAR models, and changes in user-specified models.

Clustered effects are often encountered in multiple hypothesis testing of spatial signals. In this paper, we propose a new method, termed two-dimensional spatial multiple testing (2d-SMT) procedure, to control the false discovery rate (FDR) and improve the detection power by exploiting the spatial information encoded in neighboring observations. The proposed method provides a novel perspective of utilizing spatial information by gathering signal patterns and spatial dependence into an auxiliary statistic. 2d-SMT rejects the null when a primary statistic at the location of interest and the auxiliary statistic constructed based on nearby observations are greater than their corresponding thresholds. 2d-SMT can also be combined with different variants of the weighted BH procedures to improve the detection power further. A fast step-down algorithm is developed to accelerate the search for optimal thresholds in 2d-SMT. In theory, we establish the asymptotical FDR control of 2d-SMT under weak spatial dependence. Extensive numerical experiments demonstrate that the 2d-SMT method combined with various weighted BH procedures achieves the most competitive performance in FDR and power trade-off.

Computer simulations play an important role in scientific discovery and engineering innovation. Reliable computer models enable virtual experimentation that reduces the need for costly and time-consuming physical testing. However, the credibility of such models hinges on rigorous statistical validation against real-world data. This paper develops a formal frequentist framework for both global and subdomain validation of computer models. We propose the Fourier Maximum Modulus Test (FMMT), which leverages kernel ridge regression (KRR) to estimate the discrepancy between the computer model and the physical process, followed by a frequency-domain test based on weighted generalized Fourier coefficients. The theoretical analysis establishes the asymptotic normality of these coefficients, allowing for closed-form p-values. Simulation studies and a shear-layer experiment demonstrate that FMMT achieves high power, accurate Type I error control, and strong sensitivity to localized discrepancies.

Climate data products, such as reanalysis datasets, require careful evaluation to determine their reliability. While most climate data evaluations focus on the mean and dependency of climate processes, we focus on marginal extreme behavior, including return levels that often have devastating impacts on our ecosystems and societies. In particular, we aim to identify where the two climate extreme fields exhibit different marginal behavior, by simultaneously evaluating the differences over all spatial locations through multiple testing techniques. The large variation inherited in extreme model fitting makes this evaluation more challenging than that for the mean and dependency structure. We propose a new multiple testing procedure, bivariate conditional local FDR (BiCLfdr), to efficiently detect signals from highly variable but spatially correlated hypotheses. Our method takes advantage of both the smoothness of large scale spatial variability and the local spatial correlation to enhance the power of comparing the marginal extreme distribution of two spatial extremes. We apply BiCLfdr to evaluate ERA5 reanalysis in its representation of winter precipitation extremes across the United States, using CPC-CONUS observations as reference. Our analysis identifies locations along the West Coast, the western Great Plains, and the Southeast where ERA5 appears to inadequately represent observed winter precipitation extremes.

Abstract Subcellular spatial transcriptomics resolves RNA where tissue structures extend beyond nuclei, but current aggregation strategies either impose fixed areal units or use segmented cells as anchors. This leaves process-rich and extracellular compartments difficult to measure directly. Here we report Kintsugi, a deterministic tessellation method that partitions sparse count matrices by captured-UMI density and assigns every in-tissue bin. It separates Pearson-residual gene composition from captured transcript density, preserving complete tissue representation without histology or nuclear segmentation. In mouse brain Visium HD, Kintsugi recovered nucleus-poor neuropil compartments enriched for glial and dendritic transcripts, and Xenium molecule coordinates supported dendritic mRNA localisation away from nuclei. The same representation identified nucleus-poor fibrotic scar in idiopathic pulmonary fibrosis and matrix structure in fetal cartilage. CODEX proteomics further showed that captured density was not reducible to nuclear packing alone. These results show that nucleus-poor tissue spaces contain structured, interpretable RNA biology that becomes accessible through complete tissue tessellation.

Mendelian randomization (MR) has been widely used to infer causal relationships between exposures and outcomes in epidemiological studies. However, classical MR assumptions can be violated when genetic variants are associated with outcomes through pathways other than the exposure, leading to uncorrelated and/or correlated pleiotropy. Additionally, measurement error arising from the inherent uncertainty in summary statistics obtained from large-scale genome-wide association studies can introduce bias into the causal effect estimate. To address these issues, we develop a debiased mixture inverse variance weighting ($\mathsf{dmIVW}$) method with three major advantages. First, it is capable of simultaneously handling various types of pleiotropy and eliminating the bias caused by uncertainty. Second, it can guard against distortion caused by invalid genetic variants while effectively harnessing their information. Third, our unified framework facilitates a fair comparison and combination of a series of submodels, encompassing several popular MR methods as special cases. Through real data applications, the effectiveness and robustness of $\mathsf{dmIVW}$ in estimating the causal effects of risk factors on common diseases are demonstrated.

Watermarking techniques embed statistical signals within content generated by large language models to help trace its source. Although existing methods perform well on long texts, their effectiveness significantly decreases for shorter texts. We introduce a statistical detection approach that improves the power of watermark detection, particularly in shorter texts. Our method leverages both the watermark key sequence and the next token probabilities (NTPs) to determine whether a text is generated by a large language model. We demonstrate the optimality of our approach and analyze its power properties. We also investigate an approach to estimating NTPs and extend our method to scenarios where texts face potential attacks such as substitutions, insertions, or deletions. We validate the effectiveness of our technique using texts generated by Meta-Llama-3-8B from Meta and Mistral-7B-v0.1 from Mistral AI, utilizing prompts extracted from Google's C4 dataset. In scenarios without attacks and with short text lengths, our method demonstrates approximately 65% power improvement compared to the baseline method on average.

The Bayesian Cramér-Rao bound (CRB) provides a lower bound on the mean square error of any estimator in Bayesian inference under mild regularity conditions. It benchmarks the performance of statistical estimation and can also serve as a principled metric for system design and optimization. However, it is difficult to calculate the Bayesian CRB without explicit knowledge of the prior distribution. In this paper, we introduce a novel data-driven method for Bayesian CRB estimation, leveraging state-of-the-art score estimation and deep generative modeling techniques. We show that the proposed estimator is consistent and illustrate its performance in a denoising problem.

Microbiome sequencing data are inherently sparse and compositional, with excessive zeros arising from biological absence or insufficient sampling. These zeros pose significant challenges for downstream analyses, particularly those that require log-transformation. We introduce BMDD (BiModal Dirichlet Distribution), a novel probabilistic modeling framework for accurate imputation of microbiome sequencing data. Unlike existing imputation approaches that assume unimodal abundance, BMDD captures the bimodal abundance distribution of the taxa via a mixture of Dirichlet priors. It uses variational inference and a scalable expectation-maximization algorithm for efficient imputation. Through simulations and real microbiome datasets, we demonstrate that BMDD outperforms competing methods in reconstructing true abundances and improves the performance of differential abundance analysis. Through multiple posterior samples, BMDD enables robust inference by accounting for uncertainty in zero imputation. Our method offers a principled and computationally efficient solution for analyzing high-dimensional, zero-inflated microbiome sequencing data and is broadly applicable in microbial biomarker discovery and host-microbiome interaction studies.

The Kolmogorov-Smirnov (KS) test is a widely used statistical test that assesses the conformity of a sample to a specified distribution. Its efficacy, however, diminishes with serially dependent data and when parameters within the hypothesized distribution are unknown. For independent data, parametric and nonparametric bootstrap procedures are available to adjust for estimated parameters. For serially dependent stationary data, parametric bootstrap has been developed with a working serial dependence structure. A counterpart for the nonparametric bootstrap approach, which needs a bias correction, has not been studied. Addressing this gap, our study introduces a bias correction method employing a nonparametric block bootstrap, which approximates the distribution of the KS statistic in assessing the goodness-of-fit of the marginal distribution of a stationary series, accounting for unspecified serial dependence and unspecified parameters. We assess its effectiveness through simulations, scrutinizing both its size and power. The practicality of our method is further illustrated with an examination of stock returns from the S&P 500 index, showcasing its utility in real-world applications. Supplementary materials for this article are available online.

Kolmogorov-Arnold Network (KAN) is a network structure recently proposed by Liu et al. (2024) that offers improved interpretability and a more parsimonious design in many science-oriented tasks compared to multi-layer perceptrons. This work provides a rigorous theoretical analysis of KAN by establishing generalization bounds for KAN equipped with activation functions that are either represented by linear combinations of basis functions or lying in a low-rank Reproducing Kernel Hilbert Space (RKHS). In the first case, the generalization bound accommodates various choices of basis functions in forming the activation functions in each layer of KAN and is adapted to different operator norms at each layer. For a particular choice of operator norms, the bound scales with the l_1 norm of the coefficient matrices and the Lipschitz constants for the activation functions, and it has no dependence on combinatorial parameters (e.g., number of nodes) outside of logarithmic factors. Moreover, our result does not require the boundedness assumption on the loss function and, hence, is applicable to a general class of regression-type loss functions. In the low-rank case, the generalization bound scales polynomially with the underlying ranks as well as the Lipschitz constants of the activation functions in each layer. These bounds are empirically investigated for KANs trained with stochastic gradient descent on simulated and real data sets. The numerical results demonstrate the practical relevance of these bounds.