Deformations, and diffeormophisms in particular, have played a tremendous role in the field of statistical shape analysis, as a proxy to measure and interpret differences between similar objects but with different shapes. Diffeomorphisms usually result from the integration of a flow of regular velocity fields, whose parameters have not enabled so far a full control of the local behaviour of the deformation. In this work, we propose a new mathematical and computational framework, in which diffeomorphisms are built on the combination of local deformation modules with few degrees of freedom. Deformation modules contribute to a global velocity field, and interact with it during integration so that the local modules are transported by the global diffeomorphic deformation under construction. Such modular diffeomorphisms are used to deform shapes and to provide the shape space with a sub-Riemannian metric. We then derive a method to estimate a Fréchet mean from a series of observations, and to decompose the variations in shape observed in the training samples into a set of elementary deformation modules encoding distinctive and interpretable aspects of the shape variability. We show how this approach brings new solutions to long lasting problems in the fields of computer vision and medical image analysis. For instance, the easy implementation of priors in the type of deformations offers a direct control to favor one solution over another in situations where multiple solutions may fit the observations equally well. It allows also the joint optimisation of a linear and a non-linear deformation between shapes, the linear transform simply being a particular type of modules. The proposed approach generalizes previous methods for constructing diffeomorphisms and opens up new perspectives in the field of statistical shape analysis.