In this paper, we derive error estimates of the backward Euler-Maruyama method applied to multi-valued stochastic differential equations. An important example of such an equation is a stochastic gradient flow whose associated potential is not continuously differentiable, but assumed to be convex. We show that the backward Euler-Maruyama method is well-defined and convergent of order at least $1/4$ with respect to the root-mean-square norm. Our error analysis relies on techniques for deterministic problems developed in [Nochetto, Savaré, and Verdi, Comm. Pure Appl. Math., 2000]. We verify that our setting applies to an overdamped Langevin equation with a discontinuous gradient and to a spatially semi-discrete approximation of the stochastic $p$-Laplace equation.