We calculate the decay rates for $\eta_c\to\gamma\gamma$ and $J/\psi\to\gamma \eta_c$ in lattice QCD with the effect of $u$, $d$, $s$ and $c$ quarks in the sea for the first time. Our calculations are carried out on gluon field configurations generated by the MILC Collaboration that include $2+1+1$ flavours of Highly Improved Staggered sea quarks. Valence $c$ quarks also use the Highly Improved Staggered Quark action. Extrapolation to the continuum and to physical quark masses is controlled through the use of four different lattice spacings that range from $0.015~\mathrm{fm}$ to $0.06~\mathrm{fm}$ and $u/d$ sea quarks with two masses: one-fifth that of the $s$ quark mass; and the physical $u/d$ mass. Our results are more accurate than those from previous lattice QCD calculations. This enables us to clarify considerably the theoretical picture, particularly for $\eta_c \to \gamma \gamma$ decays.