We outline a methodology for forward uncertainty quantification in systems with uncertain parameters, discontinuous model response, and a limited number of model runs. Our approach involves two stages. First we detect the discontinuity with Bayesian inference, thus obtaining a probabilistic representation of the discontinuity curve for arbitrarily distributed input parameters. Then, employing the Rosenblatt transform, we construct spectral representations of the uncertain model output, using polynomial chaos (PC) expansions on either side of the discontinuity curve, leading to an averaged PC representation of the forward model response that allows efficient uncertainty quantification. We obtain PC modes by either orthogonal projection or Bayesian inference, and argue for a hybrid approach that targets a balance between the accuracy provided by the orthogonal projection and the flexibility provided by the Bayesian inference. The uncertain model output is then computed by taking an ensemble average over PC expansions corresponding to sampled realizations of the discontinuity curve. The methodology is demonstrated on synthetic examples of discontinuous model response with adjustable sharpness and structure.