Why after transforming an object with the matrix_world of another object the reference and target object have different transformation matrices? I'm quite confused on how the matrix_world is applyed to an object.

Why the transformed object has a different transformation matrix of the reference one?
Why the parented objects do not apparently have a transformation matrix? (getting their transformation gives an error)

Here is more in detail my problem:
I'am reconstructing cerebral regions from multiple brains in a separate software (TrakEM2) and using blender (2.79) to manually align all brains togheter to a common reference space. For this I create a central plane between 3 reference points for each brain, parent all objects to that plane and align it to the reference space.
The problem is that I continually edit these brains in TrkEM2, adding details or new parts, and when importing new objects or edited ones, I need to transform them as the previoulsy imported object of the same brain. For this I use the matrix_world of one of the previously imported object as reference (anyone is good but I need to un-parent it).

the imported object goes to the right place (after an additional rotation) but the transformation matrix is different from the one I have applyed to that object. This is wierd to me, and partially a problem beacause I must remeber that these secondarily imported object should not be used as reference objects.
I tried to use the plane that I create when importing as reference, but for some reason the imported object are than shifted from where they should land (rotation ok..)

I attach a script that illustrate the firstissue, an .obj is imported, shifted and than re-imported using the first copy as a reference. As you can see, the origin and rotaion of the two copies are different, as is the printed transformation matrix.
Can anybody explain, give suggestions or send me to some documentation (for linear algebra ignorant?)
import bpy
import os

WORKING_DIR = &quot;G:\Progetto QA\RR PSD Gephryn\Export\Obj\import&quot;
file_name = &quot;RRGEF1_An_LV-pi_1226.obj&quot;
file_loc = os.path.join(WORKING_DIR, file_name)

#import object
imported_object = bpy.ops.import_scene.obj(filepath=file_loc)
obj_object = bpy.context.selected_objects[0]
bpy.context.scene.objects.active = obj_object

#Transform the object to correct default TrakEM2 orientation and scale
bpy.ops.transform.resize(value=(-0.001, 0.001, -0.001), constraint_axis=(False, False, False), constraint_orientation='GLOBAL', mirror=False, proportional='DISABLED', proportional_edit_falloff='SMOOTH', proportional_size=1)

#Translate the imported object    
bpy.ops.transform.translate(value=(2.45409, 1.31843, 2.94766), constraint_axis=(False, False, False), constraint_orientation='GLOBAL', mirror=False, proportional='DISABLED', proportional_edit_falloff='SMOOTH', proportional_size=1)

#get the transformation matrix of the transformed object
M=bpy.context.object.matrix_world
print (&quot;reference:&quot;,M)
bpy.ops.object.select_all(action='DESELECT')

#import object
imported_object = bpy.ops.import_scene.obj(filepath=file_loc)

#select imported object and apply the transformation matrix
obj_object2 = bpy.context.selected_objects[0]
bpy.context.scene.objects.active = obj_object2
me = obj_object2.data
me.transform(M)
me.update()

#Correct the transformation of the imported object
bpy.ops.transform.rotate(value=1.5708, axis=(-1, -0, -0), constraint_axis=(False, False, False), constraint_orientation='GLOBAL', mirror=False, proportional='DISABLED', proportional_edit_falloff='SMOOTH', proportional_size=1)

M2=bpy.context.object.matrix_world

print (&quot;target:&quot;,M2)
Here the updated code after following batFINGER suggestions
import bpy
import os
from mathutils import Matrix,Euler

WORKING_DIR = &quot;G:\Progetto QA\RR PSD Gephryn\Export\Obj&quot;
file_name = &quot;RRGEF1_An_LV-pi_1226.obj&quot;
file_loc = os.path.join(WORKING_DIR, file_name)

#import object
imported_object = bpy.ops.import_scene.obj(filepath=file_loc,axis_forward='Z', axis_up='-Y')
obj_object = bpy.context.selected_objects[0]
bpy.context.scene.objects.active = obj_object

obj_object.scale=(0.001,0.001,0.001)

#Translate the imported object    
bpy.ops.transform.translate(value=(2.45409, 1.31843, 2.94766), constraint_axis=(False, False, False), constraint_orientation='GLOBAL', mirror=False, proportional='DISABLED', proportional_edit_falloff='SMOOTH', proportional_size=1)

#get the transformation matrix of the transformed object
M=bpy.context.object.matrix_world

bpy.ops.object.select_all(action='DESELECT')

#re-import object
imported_object = bpy.ops.import_scene.obj(filepath=file_loc,axis_forward='Z', axis_up='-Y')

#Transform with matrix multiplication
obj_object2 = bpy.context.selected_objects[0]
obj_object2.matrix_world = M * obj_object2.matrix_world

# Adjust the transformation
obj_object2.rotation_euler=Euler((-1.5708,0,0),'XYZ')

M2=bpy.context.object.matrix_world

print (&quot;reference:&quot;,M)
print (&quot;target:&quot;,M2)