We establish a Harnack inequality for a class of quasi-linear PDE modeled on the prototype {equation*} \partial_tu= -\sum_{i=1}^{m}X_i^\ast (|\X u|^{p-2} X_i u){equation*} where $p\ge 2$, $ \ \X = (X_1,..., X_m)$ is a system of Lipschitz vector fields defined on a smooth manifold $\M$ endowed with a Borel measure $μ$, and $X_i^*$ denotes the adjoint of $X_i$ with respect to $μ$. Our estimates are derived assuming that (i) the control distance $d$ generated by $\X$ induces the same topology on $\M$; (ii) a doubling condition for the $μ$-measure of $d-$metric balls and (iii) the validity of a Poincaré inequality involving $\X$ and $μ$. Our results extend the recent work in \cite{DiBenedettoGianazzaVespri1}, \cite{K}, to a more general setting including the model cases of (1) metrics generated by Hörmander vector fields and Lebesgue measure; (2) Riemannian manifolds with non-negative Ricci curvature and Riemannian volume forms; and (3) metrics generated by non-smooth Baouendi-Grushin type vector fields and Lebesgue measure. In all cases the Harnack inequality continues to hold when the Lebesgue measure is substituted by any smooth volume form or by measures with densities corresponding to Muckenhoupt type weights.