The nonleptonic decays $\bar{B}^0 \to D^{(*)+}K^-$ and $\bar{B}^0_s \to D_s^{(*)+} \pi^-$, computed within the QCD factorization framework, are expected to be among the theoretically cleanest nonleptonic decay channels. Nevertheless, they present an intriguing tension between theoretical predictions and experimental observations. From the theory point of view this tension could be explained by underestimated theory uncertainties or beyond the Standard Model (BSM) physics. These issues have received some attention in the literature but there is no conclusive solution to the puzzle yet. In this work, we compute some of the contributions that were previously not considered. We then perform a global phenomenological study of the effective couplings, exploring the viability and possible size of a BSM explanation. The fits include constraints from the nonleptonic B-meson decay width, which we calculate at the leading order for the full set of WET operators for the first time. This study is the first one to account for simultaneous variation of up to six effective couplings. We identify two distinct modes in all fit models and discuss how future measurements can be used to break this degeneracy.