There have been various attempts to improve the reconstruction of gene regulatory networks from microarray data by the systematic integration of biological prior knowledge. Our approach follows the Bayesian paradigm where the prior knowledge is expressed in terms of energy functions, from which a prior distribution over network structures is obtained in the form of a Gibbs distribution. The hyperparameters of this distribution represent the weights associated with the prior knowledge relative to the data. We have derived and tested an MCMC scheme for sampling networks and hyperparameters simultaneously from the posterior distribution, thereby automatically learning how to trade off information from the prior and the data. We have extended this approach to a Bayesian coupling scheme for learning gene regulatory networks from a combination of related data sets that were obtained under different experimental conditions and are therefore potentially associated with different active subpathways.