This paper investigates semi-blind channel estimation for massive multiple-input multiple-output (MIMO) systems. To this end, we first estimate a subspace based on all received symbols (pilot and payload) to provide additional information for subsequent channel estimation. This additional information enhances minimum mean square error (MMSE) channel estimation. Two variants of the linear MMSE (LMMSE) estimator are formulated, where the first one solves the estimation within the subspace, and the second one uses a subspace projection as a preprocessing step. Theoretical derivations show the latter method's superior estimation performance in terms of mean square error for uncorrelated Rayleigh fading. Further, we provide asymptotic insights on how the proposed MMSE-based channel estimation strategy outperforms the unbiased Cramer-Rao bound. Subsequently, we introduce parameterizations of these semi-blind LMMSE estimators based on two different conditional Gaussian latent models, i.e., the Gaussian mixture model and the variational autoencoder. Both models learn the propagation environment's underlying channel distribution based on training data and serve as generative priors for our semi-blind channel estimation. Extensive simulations for real-world measurement data and spatial channel models show the proposed methods' superior performance compared to state-of-the-art semi-blind channel estimators in terms of MSE.