GLTF - how are joints supposed to be transformed during an animation

I'm trying to implement skeletal animation using gltf 2.0 assets.

I'm currently able to transform the skeleton and correctly render the model. The model moves as expected when a transform (for exmaple a rotation) of a joint is edited.

The problem is that as soon as I try to use the transforms from the animation sampler outputs, the skeleton is completely wrong. My testing shows that the transformation matrices of the first keyframe of the animation would match the transform of the joints in the initial pose, but they're in fact quite different ! It's not exactly clear to me where exactly these transforms are supposed to fit in the rendering algorithm.

My rendering algorithm looks roughly like this:

render_scene() {
    render_node(root_node, transform::IDENTIY)
}

render_node(node, outer_transform) {
    next_transform = outer_transform * node.transform
        
    if (node.has_skin) {
        update_joint_matrices(next_transform, node.joints)
    }
        
    if (node.has_mesh) {
        // draw calls
    }
        
    for child in node.children {
        render_node(child, next_transform)
    }
}

update_joint_matrices(outer_transform, joints) {
    world_transforms = []
    
    // Parent nodes are always processed before child nodes
    for joint in joints {
        if joint.is_root {
            world_transforms[joint] = outer_transform * joint.transform
        } else {
            world_transforms[joint] = world_transforms[joint.parent] * joint.transform
        }
    }
    
    joint_matrices = []
    
    for joint in 0..world_transforms.len() {
        joint_matrices[joint] = world_transforms[joint] * inverse_bind_matrices[joint]
    }
    
    // send joint matrices to the GPU
}

The relevant part of the vertex shader looks like this:

void main() {
    mat4 modelTransform;

    modelTransform =
        (inWeights.x * jointMatrices[int(inJoints.x)]) +
        (inWeights.y * jointMatrices[int(inJoints.y)]) +
        (inWeights.z * jointMatrices[int(inJoints.z)]) +
        (inWeights.w * jointMatrices[int(inJoints.w)]);
    }

    gl_Position = projection * view * modelTransform * vec4(inPos, 1.0);
}

Also, there's a note in the spec I don't quite understand:

Only the joint transforms are applied to the skinned mesh; the transform of the skinned mesh node MUST be ignored.

Ok, I solved the problem. The issue was that I wasn't loading quaternions correctly. Quaternions should be interpreted as raw XYZW values.