Well, I still don't see any thing wrong with the math. Even the authors rewrite the equations as sum(|A_k>|C_k>) instead of sum(|C_k>) and apply I(x)M_x instead of applying only M_x on |C_k> alone, this will not change the following equations.
I think I made a misunderstanding by the |q1>|q2> notation. Consider the entangled 2-qubits system a|00>+b|11>. The prob of the 2nd qubit is |0> is |a|^2 and the prob of 2nd qubit is |1> is |b|^2. Now apply I(x)NOT, then the system is a|01>+b|10>. The prob of 2nd qubit to be |0> is now |b|^2 and prob of 2nd qubit is |1> is |a|^2. The same can be obtained by considering the 2nd qubit only as a|0>+b|1> and apply NOT only on the second qubit, without making any operation on the first qubit and without breaking the entanglement.
I am not expert in the field, I just took a course, but I know that a mixed state is not the same as an entangled state. In their case it is entangled not mixed and so there is no need to take the partial trace. For example, consider a two subsystems entangled with each other, |q1>|q2>. I we want to apply an operator U on the second subsytem we apply I(x)U or we can simply consider the second subsystem and apply U|q2>. This is valid as long as the entanglement holds and Identity gates are assumed to be applied on the first subsystem.