The blood flow dynamics in a stenosed, subject-specific carotid bifurcation is numerically simulated using Direct Numerical Simulation (DNS) and Reynolds-averaged Navier--Stokes (RANS) equations closed with turbulence models. The former is meant to provide a term of comparison for the RANS calculations, that include classic two-equations models ($k-\epsilon$ and $k-\omega$) as well as a transitional three-equations eddy-viscosity model ($k_T-k_L-\omega$). Pulsatile inlet conditions based on in-vivo ultrasound measurements of blood velocity are used. The blood is modelled as a Newtonian fluid, and the vessel walls are rigid. The main purpose of this work is to highlight the problems related to the use of classic RANS models in the numerical simulation of such flows. The time-averaged DNS results, interpreted in view of their finite-time averaging error, are used to demonstrate the superiority of the transitional RANS model, which is found to provide results closer to DNS than those of conventional models. The transitional model is shown to possess better predictive capabilities in terms of turbulence intensity, temporal evolution of the pressure along the cardiac cycle, and the oscillatory shear index (OSI). Indeed, DNS brings to light the locally transitional or weakly turbulent state of the blood flow, which presents velocity and pressure fluctuations only in the post-stenotic region of the internal carotid artery during systole, while the flow is laminar during diastole.