We propose an efficient way to test rotational invariance in the cosmological perturbations by use of galaxy correlation functions. In symmetry-breaking cases, the galaxy power spectrum can have extra angular dependence in addition to the usual one due to the redshift-space distortion, k ·n. We confirm that, via the decomposition into not the usual Legendre basis L (k ·n) but the bipolar spherical harmonic one {Y (k)⊗Y (n)}LM , the symmetry-breaking signal can be completely distinguished from the usual isotropic one since the former yields nonvanishing L ≥ 1 modes but the latter is confined to the L = 0 one. As a demonstration, we analyze the signatures due to primordial-origin symmetry breakings such as the well-known quadrupolar-type and dipolar-type power asymmetries and find nonzero L = 2 and 1 modes, respectively. Fisher matrix forecasts of their constraints indicate that the Planck-level sensitivity could be achieved by the SDSS or BOSS-CMASS data, and an order-of-magnitude improvement is expected in a near future survey as PFS or Euclid by virtue of an increase in accessible Fourier mode. Our methodology is model-independent and hence applicable to the searches for various types of statistically anisotropic fluctuations.