The main goal of this article is to derive a two-sided estimate for hitting probabilities of a hypoelliptic stochastic differential equation (SDE) driven by fractional Brownian motion (fBM) with Hurst parameter $H\in(1/4,1)$ in terms of Newtonian-type capacities that are defined with respect to the (sub-Riemannian) control distance associated with the vector fields. As a starting point, we first establish the existence and smoothness of joint densities for the finite-dimensional distributions of the solution in the general context of hypoellitpic SDEs driven by Gaussian rough paths. We then turn to the fBM setting and derive a local upper bound for the joint density in terms of the control distance. As an application of these results, we establish our main estimate on hitting probabilities which generalises a well-known elliptic result of \cite{BNOT} to the hypoelliptic case.