The growth of large-scale spatial gene expression data requires new computational tools to extract major trends in gene expression in their native spatial context. Here, we describe an unsupervised and interpretable computational framework to (1) pre-process 3D spatial gene expression datasets by imputation of missing voxels, (2) identify principal patterns (PPs) of 3D spatial gene expression profiles using the stability-driven non-negative matrix factorization (staNMF) technique, and (3) systematically compare these PPs to known anatomical regions and ontology. This framework, referred to as osNMF (ontology discovery via staNMF), identifies PPs that are derived purely from thousands of 3D spatial gene expression profiles in the Allen Mouse Brain Atlas. These 3D PPs present stable and spatially coherent regions of the mouse brain, potentially without human labor and bias. We demonstrate that osNMF PPs offer new brain patterns that are highly correlated with combinations of expert-annotated brain regions, while also identifying a unique ontology based purely on spatial gene expression data. Compared to principal component analysis (PCA) and other clustering algorithms, our PPs exhibit better spatial coherence, more accurately match expert labeling and are more stable across multiple bootstrapped simulations. We also used osNMF to define marker genes and build putative spatial gene interaction networks. Our findings highlight the capability of osNMF to rapidly generate new atlases from a large set of spatial gene expression data without supervision and uncover novel relationships between brain regions that were difficult to discern using conventional manual approaches.