Let $(Ω,\mathcal{F},\mathbb{P})$ be a probability space and $\varphi:\ Ω\times[0,\infty)\to[0,\infty)$ be a Musielak-Orlicz function. In this article, the authors prove that the Doob maximal operator is bounded on the Musielak-Orlicz space $L^{\varphi}(Ω)$. Using this and extrapolation method, the authors then establish a Fefferman-Stein vector-valued Doob maximal inequality on $L^{\varphi}(Ω)$. As applications, the authors obtain the dual version of the Doob maximal inequality and the Stein inequality for $L^{\varphi}(Ω)$, which are new even in weighted Orlicz spaces. The authors then establish the atomic characterizations of martingale Musielak-Orlicz Hardy spaces $H_{\varphi}^s(Ω)$, $P_{\varphi}(Ω)$, $Q_{\varphi}(Ω)$, $H_{\varphi}^S(Ω)$ and $H_{\varphi}^M(Ω)$. From these atomic characterizations, the authors further deduce some martingale inequalities between different martingale Musielak-Orlicz Hardy spaces, which essentially improve the corresponding results in Orlicz space case and are also new even in weighted Orlicz spaces. By establishing the Davis decomposition on $H_{\varphi}^S(Ω)$ and $H_{\varphi}^M(Ω)$, the authors obtain the Burkholder-Davis-Gundy inequality associated with Musielak--Orlicz functions. Finally, using the previous martingale inequalities, the authors prove that the maximal Fejér operator is bounded from $H_{\varphi}[0,1)$ to $L^{\varphi}[0,1)$, which further implies some convergence results of the Fejér means; these results are new even for the weighted Hardy spaces.