A multiscalar and nonrenormalizable $B-L$ extension of the standard model (SM) with $S_4$ symmetry which successfully explains the recent observed neutrino oscillation data is proposed. The tiny neutrino masses and their hierarchies are generated via the type-I seesaw mechanism. The model reproduces the recent experiments of neutrino mixing angles and Dirac CP violating phase in which the atmospheric angle $(theta_{23})$ and the reactor angle $(theta_{13})$ get the best-fit values while the solar angle $(theta_{12})$ and Dirac CP violating phase ($delta $) belong to $3, si $ range of the best-fit value for normal hierarchy (NH). For inverted hierarchy (IH), $theta_{13}$ gets the best-fit value and $theta_{23}$ together with $de $ belongs to $1, si $ range while $theta_{12}$ belongs to $3, si $ range of the best-fit value. The effective neutrino masses are predicted to be $langle m_{ee}rangle=6.81 ,, mbox{meV}$ for NH and $langle m_{ee}rangle=48.48,, mbox{meV}$ for IH being in good agreement with the most recent experimental data.