In this paper, we propose a numerical method to solve stochastic elliptic interfaceproblems with random interfaces. Shape calculus is first employed to derive theshape-Taylor expansion in the framework of the asymptotic perturbation approach. Given themean field and the two-point correlation function of the random interface, we can thusquantify the mean field and the variance of the random solution in terms of certain ordersof the perturbation amplitude by solving a deterministic elliptic interface problem andits tensorized counterpart with respect to the reference interface. Error estimates arederived for the interface-resolved finite element approximation in both, the physical andthe stochastic dimension. In particular, a fast finite difference scheme is proposed tocompute the variance of random solutions by using a low-rank approximation based on thepivoted Cholesky decomposition. Numerical experiments are presented to validate andquantify the method.