Distorted animation of skinned model from gltf

I am exploring importing models from gltf files. I was able to render the mesh directly, and in the next iteration I am able to render a skinned model. I am creating a tree of nodes (so the ancestry is clear) and I create local_TRS matrix for each joint based on their respective translation, rotation and scale. Based on the tutorial:

jointMatrix(j) =
  globalTransformOfJointNode(j) *
  inverseBindMatrixForJoint(j);

Relevant parts of the code so far, to create joint matrices:

  applyTRS(joint, animation_matrix) {
    const parentTRS = joint.parent ? joint.parent.global_TRS : glMatrix.mat4.create();
    const animMatrix = animation_matrix[joint.jointIndex];
    glMatrix.mat4.multiply(joint.global_TRS, parentTRS, joint.local_TRS); // skin matrix
    //glMatrix.mat4.multiply(joint.global_TRS, joint.global_TRS, animMatrix); //this one ruins it

    for (const child of joint.children) {
      this.applyTRS(child, animation_matrix);
    }
    glMatrix.mat4.multiply(joint.global_TRS, joint.global_TRS, joint.InverseBindMatrix);
  }

And extraction of u_jointMat for shader uniform:

  makeJointMatrix(joint) {
    const matrix = this.skin.jointMatrix;
    const skinJoints = this.skin.skinJoints;
    const index = skinJoints.indexOf(joint.index) * 16;
    glMatrix.mat4.copy(matrix.subarray(index, index + 16), joint.global_TRS);
    for (const child of joint.children) {
      this.makeJointMatrix(child);
    }
  }

And finally, the vertex shader code:

///model_vShader///
#ifdef GL_FRAGMENT_PRECISION_HIGH
precision highp float;
#else
precision mediump float;
#endif

const int MAX_JOINTS = 64;

attribute vec4 aVertexPosition;
attribute vec3 aVertexNormal;
attribute vec2 aTextureCoord;
attribute vec4 aJoint;
attribute vec4 aWeight;

uniform mat4 uModelViewMatrix;
uniform mat4 uProjectionMatrix;
uniform mat4 uScale;
uniform mat4 uTranslate;
uniform mat4 uRotateY;
uniform mat4 u_jointMat[MAX_JOINTS];

varying vec2 vTextureCoord;
varying vec3 FragPos;
varying vec3 v_normal;

void main(void) {
    mat4 skinMat = aWeight.x * u_jointMat[int(aJoint.x)] +
        aWeight.y * u_jointMat[int(aJoint.y)] +
        aWeight.z * u_jointMat[int(aJoint.z)] +
        aWeight.w * u_jointMat[int(aJoint.w)];

    vec4 position = skinMat * aVertexPosition;
    gl_Position = uProjectionMatrix * uModelViewMatrix * uTranslate * uRotateY * uScale * position;
    vTextureCoord = aTextureCoord;
    //FragPos = vec3(aVertexPosition);
    FragPos = vec3(position); //? no effect yet?
    //v_normal = aVertexNormal;
    v_normal = vec3(skinMat * vec4(aVertexNormal, 0.0)); //? no effect yet?
}

This correctly renders the skinned model in resting pose as seen below (with ugly textured 'Cesium man'):

With the inclusion of animation data, which is interpolated and converted to corresponding TRS matrices animationMatrix:

animate(date) {
    const A = this.animations[this.animationIndex];
    const delta = ((date - this.birth) / 1000) % A.nodes[0].max;
    let keyFrameIndex = -1;
    let KFL = null;
    let keyFrameTime = 0;
    const animationMatrix = new Array(Object.keys(A.nodes).length);
    for (let nodeIndex in A.nodes) {
      const node = A.nodes[nodeIndex];
      if (!(keyFrameIndex >= 0 && KFL === node.time.length && keyFrameTime === node.time[keyFrameIndex])) {
        KFL = node.time.length;
        keyFrameIndex = binarySearchClosestLowFloat(node.time, delta, node.max / KFL);
        keyFrameTime = node.time[keyFrameIndex];
      }
      const nextKeyFrameTime = node.time[keyFrameIndex + 1];
      const timeScale = (delta - keyFrameTime) / (nextKeyFrameTime - keyFrameTime);

      const translation = glMatrix.vec3.create();
      glMatrix.vec3.lerp(translation, node.translation[keyFrameIndex], node.translation[keyFrameIndex + 1], timeScale);

      const rotation = glMatrix.quat.create();
      glMatrix.quat.slerp(rotation, node.rotation[keyFrameIndex], node.rotation[keyFrameIndex + 1], timeScale);

      const scale = glMatrix.vec3.create();
      glMatrix.vec3.lerp(scale, node.scale[keyFrameIndex], node.scale[keyFrameIndex + 1], timeScale);

      const TRS = glMatrix.mat4.create();
      glMatrix.mat4.fromRotationTranslationScale(TRS, rotation, translation, scale);
      animationMatrix[node.jointIndex] = TRS;
    }
    A.animationMatrix = animationMatrix;
    const parentJoint = this.skin.joint;
    this.applyTRS(parentJoint, A.animationMatrix);
    this.makeJointMatrix(parentJoint);
  }

If I uncomment this line in the applyTRS method:

glMatrix.mat4.multiply(joint.global_TRS, joint.global_TRS, animMatrix); //this one ruins it

The model becomes distorted, as shown below. What is not clear from the picture is that animation seems correct, feet and hands move as expected however the model I rotated unexpectedly, raised up, and relation between joins is unnatural (like the body was thrown from the building ...).

It seems like I am missing another matrix, to put it back into place. Maybe the cryptic inverse(globalTransform) from this answer:

M_joint(i) = inverse(globalTransform) * M_global(i) * inverseBindMatrix(i)

So how can I make the animation work; what am I missing?

This is the solution:

applyTRS(joint, animation_matrix) {
    const parentTRS = joint.parent ? joint.parent.global_TRS : glMatrix.mat4.create();
    const animMatrix = animation_matrix[joint.jointIndex];
    glMatrix.mat4.multiply(joint.global_TRS, parentTRS, animMatrix); // this works!
    for (const child of joint.children) {
      this.applyTRS(child, animation_matrix);
    }
    glMatrix.mat4.multiply(joint.global_TRS, joint.global_TRS, joint.InverseBindMatrix);
  }

I was wrong before, I didn't need to transform local_TRS with animation_matrix. In animation, local_TRS is TRS from animation_matrix, as seen in the code above. Now it works.