Converting bones transformation matrix from OpenGL to Blender

Question

So, I'm trying to convert some transformation matrices from OpenGL to Blender but bone poses are wrong. I assume it's transformation matrix, I'm no expert in OpenGl but it's used in this function:

glUniformMatrix4fv(glGetUniformLocation(shader_character->program, "joints"), Character::JOINT_NUM, GL_FALSE, (float*)character->joint_mesh_xform);

So I assume it's transformation matrix also I've tested it's moving bones. But when applying those values to my rig, position is weird:

Compared to OpenGL:

I've tried applying it directly with bone.matrix_basis and all axis configurations when applying rotation (like (X, -Z, Y) in Blender but none of achieved poses looks right. Also for same matrix:

<Matrix 4x4 (  0.9997,  0.0193, -0.0118, 0.0000)
            ( -0.0189,  0.9994,  0.0289, 0.0000)
            (  0.0124, -0.0287,  0.9995, 0.0000)
            (-23.8805, -7.8834, -0.4667, 1.0000)>

Blender produces different results (using matrix.to_euler()):

(-0.028720224, -0.012384685, -0.01894982)

OpenGL (using glm::extractEulerAngleXYZ()):

(0.02872, 0.012385, 0.01895)

How to achieve same pose as in OpenGL having only this transformation matrix?

Original project in OpenGL

Also linking .blend file with full data form OpenGL and rig: Link to .blend file

Edit_1:

I'm not sure but do I need to recalculate matrices earlier? I'm calculating new rotations for child bones using this code:

def forward_kinematics(self):
    for i in range(self.JOINT_NUM):
        self.joint_global_anim_xform[i] = self.joint_anim_xform[i]
        self.joint_global_rest_xform[i] = self.joint_rest_xform[i]
        j = int(self.joint_parents[i])
        while j != -1:
            self.joint_global_anim_xform[i] = self.joint_global_anim_xform[i] @ self.joint_anim_xform[j]
            self.joint_global_rest_xform[i] = self.joint_global_rest_xform[i] @ self.joint_rest_xform[j]
            j = int(self.joint_parents[j])
​
        self.joint_mesh_xform[i] = self.joint_global_rest_xform[i].inverted() @ self.joint_global_anim_xform[i]

And I'm applying rotations using this code:

def apply_rotations(self, frame :int):
    for bone in self.armature.pose.bones:
        ind = self.bones.get(bone.name, -1)
        if ind != -1:
            bone.matrix_basis = mat.Matrix(self.joint_mesh_xform[ind])

Matrices used here have exact same values as ones used in OpenGL version

Edit_2:

Class definition for ealier code:

class Character:
    def __init__(self):
        self.JOINT_NUM = 31
        self.joint_anim_xform = np.empty(self.JOINT_NUM, dtype=mat.Matrix)
        self.joint_rest_xform = np.empty(self.JOINT_NUM, dtype=mat.Matrix)
        self.joint_mesh_xform = np.empty(self.JOINT_NUM, dtype=mat.Matrix)
        self.joint_global_rest_xform = np.empty(self.JOINT_NUM, dtype=mat.Matrix)
        self.joint_global_anim_xform = np.empty(self.JOINT_NUM, dtype=mat.Matrix)
        self.joint_parents = [0] * self.JOINT_NUM

rest_xform are default bones positions, anim_xform are calculated bones positions

Thanks in advance for help!