  We study general parabolic equations of the form $u_t = div A(x,t, u,D u) +
div(|F|^{p-2} F)+ f$ whose principal part depends on the solution itself. The
vector field $A$ is assumed to have small mean oscillation in $x$, measurable
in $t$, Lipschitz continuous in $u$, and its growth in $Du$ is like the
$p$-Laplace operator. We establish interior Calder\'on-Zygmund estimates for
locally bounded weak solutions to the equations when $p>2n/(n+2)$. This is
achieved by employing a perturbation method together with developing a
two-parameter technique and a new compactness argument. We also make crucial
use of the intrinsic geometry method by DiBenedetto \cite{D2} and the maximal
function free approach by Acerbi and Mingione \cite{AM}.
