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, $ \hat{k} \cdot \hat{n}$. We confirm that, via the decomposition into not the usual Legendre basis ${\cal L}_\ell(\hat{k} \cdot \hat{n})$ but the bipolar spherical harmonic one $\{Y_{\ell}(\hat{k}) \otimes Y_{\ell'}(\hat{n})\}_{LM}$, the symmetry-breaking signal can be completely distinguished from the usual isotropic one since the former yields nonvanishing $L \geq 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.