Bayesian networks are graphical models to represent the deterministic and probabilistic relationships between variables within the Bayesian framework. The knowledge of all variables can be updated using new information about some of the variables. The Bayesian Generalized Linear Least Squares method can be regarded as an inference method for Bayesian networks of variables with multivariate normal priors and linear relationships between them. We show that relying explicitly on the Bayesian network interpretation enables large scale inference and gives more flexibility in incorporating prior assumptions and constraints into the nuclear data evaluation process, such as the constraints that some cross sections equal linear combinations of other cross sections and that all cross sections must be non-negative. The latter constraint is accounted for by a nonlinear transformation and therefore we also discuss inference in Bayesian networks with non-linear relationships between variables. Using Bayesian networks, the evaluation process yields more detailed information, such as posterior estimates and uncertainties of all statistical and systematic errors associated with the experiments. We further elaborate on a sparse Gaussian process construction that can be well integrated into the Bayesian network framework and applied to, e.g., the modeling of energy-dependent model parameters, model deficiencies of the physics model or energy-dependent systematic errors of experiments. We present three proof-of-concept examples that emerged in the context of the neutron data standards project and in the ongoing international evaluation efforts of 56 Fe. In the first example we demonstrate the modelization and explicit estimation of relative energy-dependent error components associated with experimental datasets. Then we show that Bayesian networks in combination with the outlined Gaussian process construction may be applied to an evaluation of 56 Fe in the energy range between one and two MeV, where it is difficult to obtain satisfactory evaluations by R-Matrix and nuclear model fits. Finally, we present a model-based evaluation of 56 Fe between 5 MeV and 30 MeV with a consistent and statistically sound treatment of model deficiencies. The R scripts to reproduce the Bayesian network examples and the nucdataBaynet package for Bayesian network modeling and inference have been made publicly available.