In the absence of symmetry assumptions most numerical relativity simulations adopt Cartesian coordinates. While Cartesian coordinates have some desirable properties, spherical polar coordinates appear better suited for certain applications, including gravitational collapse and supernova simulations. Development of numerical relativity codes in spherical polar coordinates has been hampered by the need to handle the coordinate singularities at the origin and on the axis, for example by careful regularization of the appropriate variables. Assuming spherical symmetry and adopting a covariant version of the BSSN equations, Montero and Cordero-Carri\'on recently demonstrated that such a regularization is not necessary when a partially implicit Runge-Kutta (PIRK) method is used for the time evolution of the gravitational fields. Here we report on an implementation of the BSSN equations in spherical polar coordinates without any symmetry assumptions. Using a PIRK method we obtain stable simulations in three spatial dimensions without the need to regularize the origin or the axis. We perform and discuss a number of tests to assess the stability, accuracy and convergence of the code, namely weak gravitational waves, "hydro-without-hydro" evolutions of spherical and rotating relativistic stars in equilibrium, and single black holes.