The author considers the \(p\)-harmonic system \[ \text{div} \left( \langle A(x)Du(x),Du(x)\rangle^{(p-2)/2} A(x)Du(x) \right) = \text{div}\left( \sqrt{A(x)}F(x) \right) \tag \(*\) \] where \(A(x)\) is a positive definite matrix whose entries have bounded mean oscillation (BMO) and \(F\) is a given matrix field. The main result is the following a priori estimate. If \(|r-2|0\) depending on the BMO norm of \(\sqrt A\), then there is a \(\delta>0\) such that \[ \|\sqrt{A} Du \|_r^r \leq C \|F \|_{r/(p-1)}^{r/(p-1)} \] whenever \(|p-r|<\delta\) and \(u\) is a very weak solution of \((*)\). The proof uses the existence and uniqueness results of the author [Potential Analysis 15, No. 3, 285-299 (2001)] for linear systems with BMO coefficients combined with nonlinear commutator relations.