This paper presents a comprehensive methodology that combines uncertainty quantification, uncertainty propagation, and design optimization using a Bayesian framework. The epistemic uncertainty due to input data uncertainty is considered. Two types of uncertainty models for input variables and/or their distribution parameters are addressed: 1) uncertainty modeled as family of distributions, and 2) uncertainty modeled as interval data. A Bayesian approach is adopted to update the uncertainty models, where the likelihood functions are constructed using limited experimental data. Global sensitivity analysis, which previously only considered aleatory inputs in the context of probabilistic representation, is extended in this paper to quantify the contributions of both aleatory and epistemic uncertainty sources for multioutput problems using an auxiliary variable approach. Gaussian process surrogate modeling is employed to replace the expensive physics models and improve the computational efficiency. A previously developed bias-minimization technique, which only dealt with single-output functions, is extended to reduce the surrogate model error for a multioutput function. A decoupled robustness-based design optimization framework is developed to include both aleatory and epistemic uncertainties. The proposed methodology is illustrated using the NASA Langley Research Center’s multidisciplinary uncertainty quantification challenge problem.