This paper proposes a Bayesian nonparametric framework for clustering axially symmetric data. Our approach is based on a Dirichlet processes mixture model with Watson distributions, which can also be considered as the infinite Watson mixture model. In this paper, first, we extend the finite Watson mixture model into its infinite counterpart based on the framework of truncated Dirichlet process mixture model with a stick-breaking representation. Second, we propose a coordinate ascent mean-field variational inference algorithm that can effectively learn the parameters of our model with closed-form solutions; Third, to cope with a massive data set, we develop a stochastic variational inference algorithm to learn the proposed model through the method of stochastic gradient ascent; Finally, the proposed nonparametric Bayesian model is evaluated through simulated axially symmetric data sets and a real-world application, namely, gene expression data clustering.