We create data-driven maps of transcriptomic anatomy with a probabilistic framework for unsupervised pattern discovery in spatial gene expression data. Convolved negative binomial regression is used to find patterns which correspond to cell types, microenvironments, or tissue components, and that consist of gene expression profiles and spatial activity maps. Expression profiles quantify how strongly each gene is expressed in a given pattern, and spatial activity maps reflect where in space each pattern is active. Arbitrary covariates and prior hierarchies can be specified to support and leverage complex experimental designs.Transcriptomic patterns in mouse brain and olfactory bulb correspond to neuroanatomically distinct cell layers. Moreover, batch effects are successfully addressed, leading to consistent pattern inference for multi-sample analyses. On this basis, we identify known and uncharacterized genes that are spatially differentially expressed in the hippocampal field between Ammon's horn and the dentate gyrus.The analysis of spatially stratified, transcriptome-wide data [1][2][3] poses additional challenges compared to classical analysis of bulk RNA sequencing (RNA-Seq) samples. In the classical setting, annotated covariates dictate how samples are to be grouped and compared. These covariates are typically known, and often controlled for, as is the case when performing differential gene expression (DGE) analysis for bulk sequencing count data in DESeq2 [4]. In contrast, the covariates determining gene expression in space are often unknown and can change both gradually or abruptly. For example, the number of infiltrating immune cells per unit area often If the cell types underlying a sample are well characterized, it becomes analytically possible to determine mixing proportions of the cell types' known expression profiles in spatial gene expression data. But in the general case, there exists the dual discovery problems of expression profiles and spatial distributions.Two recently introduced statistical tests, Trendsceek and SpatialDE, quantify the extent of spatial variation of individual genes' expression [5,6]. Transcending individual genes, the "automatic expression histology" feature of SpatialDE implements a gene clustering approach to hidden pattern discovery. Such clustering formulations, however, appear challenged by genes that participate in multiple, spatially overlapping expression programs, and grade-of-membership formulations seem more appropriate for hidden pattern discovery.Joint analysis of multiple data sets substantially increases statistical power and sensitivity but is challenging due to the presence of batch effects. Both nonparametric [7,8] and model-based approaches are used to address such effects in the analysis of single cell RNA sequencing (scRNA-Seq) data. The model-based methods perform regression for the count data based on known sample-level covariates and allow for the discovery of unknown ones.Log-normal expression models are used for bulk [9,10] and single cell [11,12] R...