Summary: We propose a novel methodology based on a Bayesian Gaussian conditional random field model for elegantly learning the conditional dependence structures among multiple outcomes, and between the outcomes and a set of covariates simultaneously. Our approach is based on a Bayesian hierarchical model using a spike and slab Lasso prior. We investigate the corresponding maximum a posteriori (MAP) estimator that requires dealing with a nonconvex optimization problem. In spite of the nonconvexity, we establish the statistical accuracy for all points in the high posterior region, including the MAP estimator, and propose an efficient EM algorithm for the computation. Using simulation studies and a real application, we demonstrate the competitive performance of our method for the purpose of learning the dependence structure.