Recent advances in spatial transcriptomics (ST) and cost reductions have enabled large-scale multi-slice ST data generation, enhancing the statistical power to detect subtle biological signals. However, cross-slice inconsistencies and data quality variability present significant analytical challenges. To overcome these limitations, we developed MUSE, a computational framework designed for multislice joint embedding, spatial domain identification, and gene expression imputation. Specifically, MUSE integrates a two-module architecture to ensure robust cross-slice alignment and data harmonization. The alignment module models each slice as a graph and employs optimal transport to align cells across slices while preserving spatial continuity. The optimization module further refines integration by incorporating an alignment loss, allowing lower-quality data to leverage structural information from higher-quality slices. Additionally, MUSE generates virtual neighbors from aligned cells, enriching contextual information and mitigating data sparsity. These design principles enable seamless integration with existing single-slice methods, extending their applicability to multi-slice ST analysis. To comprehensively evaluate its performance, we applied MUSE to 12 real and 48 simulated datasets spanning a range of data qualities. Across all metrics, MUSE consistently outperformed existing methods in cross-slice consistency, spatial domain identification, and gene expression imputation. To promote accessibility and adoption, we provide MUSE as an open-source software package. As multi-slice ST datasets become increasingly prevalent, MUSE provides a robust and extensible framework designed to effectively integrate growing numbers of slices, thereby advancing the analysis of tissue architectures and spatial gene expression in complex biological systems.