Recent discovery of kilohertz quasi-periodic brightness oscillations of low mass X-ray binaries (LMXBs) has attracted attention to highly relativistic periodic motion near accreting neutron stars. Most models proposed so far involve (almost) inertial motion in the vicinity of the stars' innermost stable circular orbits. In the present paper we study general-relativistic circular and eccentric orbits around spinning neutron stars assuming the orbits are slightly tilted with respect to the stars' equatorial planes. We develop analytical and numerical techniques for integrating bound timelike geodesics in fully relativistic neutron star spacetimes obtained by modern numerical codes. We use equations of state of neutron star matter that span a broad range of stiffness, while the explored range of masses M > 1.7 M(solar) and spin frequencies nu_s < 600 Hz is motivated by the observations of LMXBs. We investigate the general-relativistic effects of periastron advance and nodal precession in the strong gravitational fields of rotating neutron stars and compare quantitatively the associated orbital frequencies with the more readily obtainable frequencies of orbits around Kerr black holes on the one hand, and low-order post-Newtonian (PN) expansions, on the other. While Kerr results approximate the periastron advance frequency much better than the PN expressions, the retrograde torque caused by the high quadrupole moments of the rotating stars clearly favours the PN approximation in the case of nodal precession. The methods developed in the present paper are used in the companion paper to test recent hypotheses regarding the origin of quasi-periodic oscillations in accreting neutron star sources.