We present and analyze several numerical methods for the discretization of the nonlinear Dirac equation in the nonrelativistic limit regime, involving a small dimensionless parameter \(0<\varepsilon\ll 1\) which is inversely proportional to the speed of light. In this limit regime, the solution is highly oscillatory in time, i.e. there are propagating waves with wavelength \(O(\varepsilon^2)\) and \(O(1)\) in time and space, respectively. We begin with four frequently used finite difference time domain (FDTD) methods and establish rigorously error estimates for the FDTD methods, which depend explicitly on the mesh size \(h\) and time step \(\tau\) as well as the small parameter \(0<\varepsilon\le 1\). Based on the error bounds, in order to obtain `correct' numerical solutions in the nonrelativistic limit regime, i.e. \(0<\varepsilon\ll 1\), the FDTD methods share the same \(\varepsilon\)-scalability: \(\tau=O(\varepsilon^3)\) and \(h=O(\sqrt{\varepsilon})\). Then we propose and analyze two numerical methods for the discretization of the nonlinear Dirac equation by using the Fourier spectral discretization for spatial derivatives combined with the exponential wave integrator and time-splitting technique for temporal derivatives, respectively. Rigorous error bounds for the two numerical methods show that their \(\varepsilon\)-scalability is improved to \(\tau=O(\varepsilon^2)\) and \(h=O(1)\) when \(0<\varepsilon\ll 1\) compared with the FDTD methods. Extensive numerical results are reported to confirm our error estimates.