Phoretic motion of particle clusters suspended in a liquid or a gas is observed when an external gradient of some physical property (temperature, species concentration, electrostatic potential, etc.) is applied to the system. In this study, a theory of phoretic motion of clusters of arbitrary number of spherical particles is built for example of thermophoresis in the near-continuum regime. Using a multipole expansion of the temperature and flow velocity in a series of spherical harmonics and formulating the Lamb's boundary conditions for velocity on the cluster surface, the problem is reduced to solution of the infinite system of linear equations for the expansion coefficients. For practical applications, the system is truncated and the truncation level is the only parameter, determining the accuracy of obtained numerical solution. The method is characterized by good convergence, as shown by comparing predictions of the described theory with available theoretical results. Based on the developed theory, influence of shape and size of various particle clusters on their thermophoretic velocity is established. In another application, the electrophoretic behavior of binary clusters, consisting of spherical colloidal particles with different signs of zeta-potentials, is investigated. It is found that ensemble-averaged electrophoretic velocity in general case can be adequately described using the Smoluchowski model with average value of zeta potential in the cluster.