In this work, we studied methanol and its deuterated derivatives (CH3OH, CH3OD, CD3OH, CD3OD) by NIR spectroscopy and anharmonic quantum chemical calculations. Vibrational bands corresponding to up to three quanta transitions (first and second overtones, binary and ternary combination modes) were predicted by the use of the VPT2 route. The accuracy of prediction of NIR modes was evaluated through density functional theory (DFT) with selected density functionals and basis sets. On the basis of the theoretical NIR spectra, detailed band assignments for all studied molecules were proposed. It was found that the pattern of bands in NIR spectra of deuterated methanols can be used for identification of isotopically equalized forms. Calculations of NIR spectra of all possible forms of CXXXOX (X = H, D) molecules demonstrated that the isotopic contamination can be identified due to a coexistence of bands specific to OH and OD groups. Also, bands from partially deuterated methyl groups can be distinguished in NIR spectra. Since the VPT2 framework is known to be sensitive to inaccuracy in the case of highly anharmonic modes, we obtained an independent insight by numerical solving of the time-independent Schrödinger equation corresponding to the O-X stretching mode scanned within -0.4 to 2.0 Å over a dense grid of 0.005 Å. This way the energies of vibrational levels of the CX1X2X3OX4 (X = H, D) isotopomers and the corresponding transition frequencies were obtained with high accuracy (<0.1 cm-1). The change in normal coordinate influences the reduced mass of the oscillator and thus its frequency. Our results lead to a conclusion that the effect of deuterization of the methyl group introduces a very specific and consistent frequency shift of the first overtone of the O-X stretching mode depending on the substitution of X1, X2, or X3 positions (<2 cm-1). However, the pattern of this shift is not reproduced accurately and is also largely overestimated by VPT2 calculations.