Here are some small facts about Quaternionic Geometry. First let us leave metrics aside, for the moment. There are 2 analogues of complex manifolds in the quaternionic world, namely triholomorphic manifolds and quaternionic manifolds. They are very different. Triholomorphic manifolds admit a 2-sphere of complex structures. When I write complex structure, I mean integrable almost complex structure. On the other hand, quaternionic manifolds are manifolds admitting an almost quaternionic structure satisfying an integrability condition. They may not admit a 2-sphere of almost complex structures globally, though they do locally. But these almost complex structures need not be integrable.
As an example of a triholomorphic manifold, think of $\mathbb{H}^n$ for instance. An example of a quaternionic manifold would be $\mathbb{H}P^n$, which does not even admit a globally defined almost complex structure.
The metric analogues of these 2 classes of manifolds, would be hyperkähler manifolds and quaternion-kähler manifolds. Their definitions can be found in many places (wikipedia for instance, or Besse's Einstein manifolds, or books or lecture notes by say, Salamon for instance).
I hope this will point you in the right direction. Note that hyperkähler manifolds have holonomy which is a subgroup of $Sp(k)$, where the dimension of the manifold is $4k$, while the holonomy of a quaternion-kähler manifold of dimension $4k$ is a subgroup of $Sp(k)\times_{\mathbb{Z}_2} Sp(1)$.
Edit: I will answer your question directly, for hyperkähler manifolds, where the $J_u$ are integrable and globally defined. If $u,v,w$ is a cyclic permutation of $1,2,3$, then:
$g_{st}(J_u)^s_{\,i}(J_v)^t_{\,j} = -g_{st}(J_u)^s_{\,i}(J_u)^t_{\,k}(J_w)^k_{\,j} = -g_{ik}(J_w)^k_{\,j}$
So the general answer is:
$g_{st}(J_u)^s_{\,i}(J_v)^t_{\,j} = g_{ij} \delta_{uv} - \epsilon_{uvw}g_{ik}(J_w)^k_{\,j}$
where summation over repeated indices was implicitly assumed.
Edit: the OP asked me in the comment about the case where $J_u$ is denoted by $(J_u)_i^{\,s}$ say (the "transpose" convention, so to speak). We repeat the first calculation here, with this alternate convention.
If $u,v,w$ is a cyclic permutation of $1,2,3$, then:
$g_{st}(J_u)^{\,s}_i(J_v)^{\,t}_j = g_{st}(J_u)^{\,s}_i(J_w)^{\,k}_j(J_u)^{\,t}_k = g_{ik}(J_w)^{\,k}_j$
So the general answer with this alternate convention is:
$g_{st}(J_u)^{\,s}_i(J_v)^{\,t}_j = g_{ij} \delta_{uv} + \epsilon_{uvw}g_{ik}(J_w)^{\,k}_j$