Optimization models involving quaternion matrices are widely used in color image process and other engineering areas. These models optimize real functions of quaternion matrix variables. In particular, 0-norms and rank functions of quaternion matrices are not continuous at all. Yet calculus with derivatives, subdifferentials and generalized subdifferentials of such real functions is needed to handle such models. In this paper, we introduce first and second order derivatives and establish their calculation rules for such real functions. Our approach is consistent with the subgradient concept for norms of quaternion matrix variables, recently introduced in the literature. We develop the concepts of generalized subdifferentials of lower semi-continuous functions of quaternion matrices, and use them to analyze the optimality conditions of a sparse low rank color image denoising model. We introduce R-product for two quaternion matrix vectors, as a key tool for our calculus. We show that the real representation set of low-rank quaternion matrices is closed and semi-algebraic. We also establish first order and second order optimality conditions for constrained optimization problems of real functions in quaternion matrix variables.