Firstly, I will answer a question similar to yours. The essence of quantum field theory is that we typically consider small perturbations around free theories. That is, given a set of $N_b$ Bosonic fields $\phi_i$ and $N_f$ Fermionic fields $\psi_i$ (let us consider only spin $0$ and spin $1/2$), then the typical action considered is given by
$$S=\sum_{i=1}^{N_b}S_0[\phi_i]+\sum_{i=1}^{N_f}S_0[\psi_i]+\sum_{n=1}^{\infty}\epsilon^nS^{(n)}_{\text{int}}[\phi,\psi]$$
Where $\epsilon$ is a small parameter around which we perturb. The expansion of the fields in terms of fourier modes is only truly valid at $\epsilon=0$, at which point there is no interaction among the fields, and we the full action is simply the sum of several free actions of different particles. Thus, each is quantized independently since the time-evolution operator has no chance of mixing the fields.
Thus, when $\epsilon=0$, there is no way that the Fourier modes (creation and annihilation operators) of different fields could (anti)commute to give a nonzero result. (I should be careful saying this. Let us define a bracket $[.,.\}$ such that it yields a commutator for two Bosonic operators, an anticommutator for two Fermionic operators, and a commutator for one Bosonic and Fermionic operator. Then this statement is that, given two Fourier modes $a$ and $b$ from different fields, we should have $[a,b\}=0$.) Then, interactions are given by the fact that the interaction action $S_{\text{int}}$ contains products of Fourier modes from different fields.
Now, to answer the first part of your question, about the commutation relations of complex scalar fields. Let us consider the free particle action
$$S_0=\partial_{\mu}\varphi^{*}\partial^{\mu}\varphi$$
This can be reduced to the sum of scalar particle actions by introducing $\varphi_1=(\varphi+\varphi^{\dagger})/\sqrt{2}$ and $\varphi_2=(\varphi-\varphi^{\dagger})/\sqrt{2}\,i$, which means that $\varphi_1$ and $\varphi_2$ are real. The action here becomes
$$S_0=\frac{1}{2}\partial_{\mu}\varphi_1\partial^{\mu}\varphi_1+\frac{1}{2}\partial_{\mu}\varphi_2\partial^{\mu}\varphi_2$$
Given the decomposition you gave, we can introduce $c_{\textbf{k}}\equiv(a_{\textbf{k}}+b_{\textbf{k}})/\sqrt{2}$ and $d_{\textbf{k}}\equiv(a_{\textbf{k}}-b_{\textbf{k}})/\sqrt{2}\,i$, then we can expand our fields as
$$\varphi_1(x)=\int\frac{\mathrm{d}^D\textbf{k}}{(2\pi)^D}\frac{1}{\sqrt{2\omega_{\textbf{k}}}}\left[c_{\textbf{k}}e^{ik\cdot x}+c^{\dagger}_{\textbf{k}}e^{-ik\cdot x}\right]$$
$$\varphi_2(x)=\int\frac{\mathrm{d}^D\textbf{k}}{(2\pi)^D}\frac{1}{\sqrt{2\omega_{\textbf{k}}}}\left[d_{\textbf{k}}e^{ik\cdot x}+d^{\dagger}_{\textbf{k}}e^{-ik\cdot x}\right]$$
Now, because of the analysis from earlier, we know that the $c$ and $d$ operators commute. Furthermore the canonical quantization gives us
$$\left[c_{\textbf{k}},c^{\dagger}_{\textbf{p}}\right]=\left[d_{\textbf{k}},d^{\dagger}_{\textbf{p}}\right]=(2\pi)^D\delta(\textbf{k}-\textbf{p})$$
With all other relations vanishing. From this, it is not difficult to show that the $a$ and $b$ operators commute with each other.
There is a matter when this is much more subtle, and the quantization by introducing Fourier modes is not useful. If, in the action at the top of this answer does not have $\epsilon\ll 1$, then we can no longer think of $S$ as a deformation of a free quantum field theory. The classical equations of motion will no longer be those of a free particle, and thus, since the interactions will depend on combinations of different fields, we will no longer have the translational invariance that made the Klein-Gordon equation soluble by Fourier decomposition.
Now, I will show an example of when a strongly-interacting field theory exists in nature (QCD). Consider the action for $SU(N)$ Yang-Mills theory (if you haven't studied non-abelian gauge theory yet, which seems to be the case, the analysis is still accessible), given by
$$S=\frac{1}{4g^2}\int\mathrm{d}^{D+1}x\,\text{Tr}\left[F^{\mu\nu}F_{\mu\nu}\right]$$
Where $F_{\mu\nu}$ is a tensor of $N\times N$ matrices given by
$$F_{\mu\nu}=\partial_{\mu}A_{\nu}-\partial_{\nu}A_{\mu}+[A_{\mu},A_{\nu}]$$
(Note that we are working in the Euclidean signature, where the metric is given by $\delta_{\mu\nu}$). The $A_{\mu}$ are vector fields (traditionally gluons) given by $N\times N$ matrices living in the Lie Algebra of $SU(N)$. It is clear that this action has cubic and quartic interactions among the gluons. Now, if we make the field redefinitions $A\to gA$, then cubic terms in the action are weighted by $g$ and the quartic terms are weighted by $g^2$. Thus, everything is fine and dandy if $g\ll 1$. However, if $g$ is large, then the interaction terms can no longer be considered perturbations, and Fourier decomposition is no longer useful. This is exactly the case in QCD, where the coupling constant is large at low energies.
I hope this helped!
