We present a general framework for the rigorous numerical analysis of time-fractional nonlinear parabolic partial differential equations, with a fractional derivative of order $\alpha\in(0,1)$ in time. It relies on three technical tools: a fractional version of the discrete Grönwall type inequality, discrete maximal regularity, and regularity theory of nonlinear equations. We establish a general criterion for showing the fractional discrete Grönwall inequality and verify it for the L1 scheme and convolution quadrature generated by backward difference formulas. Further, we provide a complete solution theory, e.g., existence, uniqueness, and regularity, for a time-fractional diffusion equation with a Lipschitz nonlinear source term. Together with the known results of discrete maximal regularity, we derive pointwise $L^2(\Omega)$ norm error estimates for semidiscrete Galerkin finite element solutions and fully discrete solutions, which are of order $O(h^2)$ (up to a logarithmic factor) and $O(\tau^\alpha)$, respectively, without any extra regularity assumption on the solution or compatibility condition on the problem data. The sharpness of the convergence rates is supported by the numerical experiments.