With the growing amount and diversity of intermediate omics data complementary to genomics (e.g., DNA methylation, gene expression, and protein abundance), there is a need to develop methods to incorporate intermediate omics data into conventional genomic evaluation. The omics data helps decode the multiple layers of regulation from genotypes to phenotypes, thus forms a connected multi-layer network naturally. We developed a new method named NN-LMM to model the multiple layers of regulation from genotypes to intermediate omics features, then to phenotypes, by extending conventional linear mixed models (“LMM”) to multi-layer artificial neural networks (“NN”). NN-LMM incorporates intermediate omics features by adding middle layers between genotypes and phenotypes. Linear mixed models (e.g., pedigree-based BLUP, GBLUP, Bayesian Alphabet, single-step GBLUP, or single-step Bayesian Alphabet) can be used to sample marker effects or genetic values on intermediate omics features, and activation functions in neural networks are used to capture the nonlinear relationships between intermediate omics features and phenotypes. NN-LMM had significantly better prediction performance than the recently proposed single-step approach for genomic prediction with intermediate omics data. Compared to the single-step approach, NN-LMM can handle various patterns of missing omics measures, and allows nonlinear relationships between intermediate omics features and phenotypes. NN-LMM has been implemented in an open-source package called “JWAS”.