Aligning an object to the surface of a sphere while maintaining forward direction

Using OpenGL as a rendering engine I'm attempting to generate a rotation matrix that will sit a mesh on the surface of a sphere at a given position. Initially I'd attempted to do this by calculating the spheres surface normal at the given position and using that as my meshes up vector.

This worked, but things are complicated by the fact that I'd also like to be able to rotate the mesh with absolute yaw, pitch and roll values - where a yaw of 0 will always leave the mesh pointing upwards (0, 1, 0) in a northern direction.

I have a mental block on how to achieve this and would appreciate any input.

Update with solution provided by Ripi2:

let earthAxisOfRotation = Vector3(0, 1, 0).normalized() // south -> north pole is `up` in world space
let surfaceNormal = position.normalized() // surface normal of the given position on sphere surface

let east = earthAxisOfRotation.cross(surfaceNormal).normalized()
let north = surfaceNormal.cross(east)

let rotationMatrix = Matrix4(
    east.x, east.y, east.z, 0,
    surfaceNormal.x, surfaceNormal.y, surfaceNormal.z, 0,
    -north.x, -north.y, -north.z, 0,
    0, 0, 0, 1
)

Solution

Let's work on an "Earth". Its axis of rotation (r) goes from south pole to north pole. A normal (n) to the surface goes "out" of Earth.

The cross product e= r x n gives a vector heading east.
The cross product d= n x e gives a vector heading north.

With normalized versions of e, n, d and translating to the point in the surface you can construct a lookAt matrix (see here).
up = n and f= center-eye = d and s = e. Notice you may want to change f sign because this is not a real 'lookAt' transform, it doesn't look at, but just locates.

Now, if your mesh can be considered as another sphere concentric with the first sphere, all rotations must be over the shared center of both spheres. The axis of rotation, keeping always the "heading north" can be the "Earth axis" (r) and the "heading east" (e). When using 'e' you lose somehow north, vector goes "inside" the sphere. Angles can be computed from cross product, due to if c= a x b then |c| = |a|·|b|·sin(beta).