A recursion procedure for the analytical generation of hyperspherical harmonics for tetraatomic systems, in terms of row-orthonormal hyperspherical coordinates, is presented. Using this approach and an algebraic Mathematica program, these harmonics were obtained for values of the hyperangular momentum quantum number up to 30 (about 43.8 million of them). Their properties are presented and discussed. Since they are regular at the poles of the tetraatomic kinetic energy operator, are complete, and are not highly oscillatory, they constitute an excellent basis set for performing a partial wave expansion of the wave function of the corresponding Schrodinger equation in the strong interaction region of nuclear configuration space. This basis set is, in addition, numerically very efficient and should permit benchmark-quality calculations of state-to-state differential and integral cross sections for those systems.