We present the first effective-one-body (EOB) model for generic-spins quasi-circular black-hole--neutron-star (BHNS) inspiral-merger-ringdown gravitational waveforms (GWs). Our model is based on a new numerical-relativity (NR) informed expression of the BH remnant and its ringdown. It reproduces the NR $(\ell,m)=(2,2)$ waveform with typical phase agreement of ${\lesssim0.5}\,$rad (${\lesssim1}\,$rad) to merger (ringdown). The maximum (minimum) mismatch between the $(2,2)$ and the NR data is 4% (0.6%). Higher modes (HMs) $(2,1)$, $(3,2)$, $(3,3)$, $(4,4)$ and $(5,5)$ are included and their mismatch with the available NR waveforms are up to (down to) a 60% (1%) depending on the inclination. Phase comparison with a 16 orbit precessing simulation shows differences within the NR uncertainties. We demonstrate the applicability of the model in GW parameter estimation by perfoming the first BHNS Bayesian analysis with HMs (and non-precessing spins) of the event GW190814, together with new $(2,2)$-mode analysis of GW200105 and GW200115. For the GW190814 study, the inclusion of HMs gives tighter parameter posteriors. The Bayes factors of our analyses on this event show decisive evidence for the presence of HMs, but no clear preference for a BHNS or a binary black hole (BBH) source. Similarly, we confirm GW200105 and GW200115 show no evidence for tidal effects.