We develop a method to compute low-eccentricity initial data of black hole--neutron star binaries in the puncture framework extending previous work on other types of compact binaries. In addition to adjusting the orbital angular velocity of the binary, the approaching velocity of a neutron star is incorporated by modifying the helical Killing vector used to derive equations of the hydrostationary equilibrium. The approaching velocity of the black hole is then induced by requiring the vanishing of the total linear momentum of the system, differently from the case of binary black holes in the puncture framework where the linear momentum of each black hole is specified explicitly. We successfully reduce the orbital eccentricity to <~0.001 by modifying the parameters iteratively using simulations of ~3 orbits both for nonprecessing and precessing configurations. We find that empirical formulas for binary black holes derived in the excision framework do not reduce the orbital eccentricity to ~0.001 for black hole--neutron star binaries in the puncture framework, although they work for binary neutron stars.