The QCD factorization (QCDF) method usually contains infrared divergences which introduce large model dependence to its predictions on charmless B decays. The amplitudes of charmless B decays can be decomposed into "tree" and "penguin" parts which are conventionally defined, not from the topology of the dominant diagrams, but through their associated CKM factors V * ub V uq and V * tb V tq , respectively, with q = d, s. We find that for B d,s → π ∓ K ± decays, the "tree" amplitude can be well estimated in QCDF with small errors, as the endpoint singularities have been canceled to a large extent. With this as the only input from QCDF and combined with flavor symmetry, the branching ratio of B s → π + K − is estimated to be significantly larger than the CDF measurement.This contradiction could be solved if the form factor F BsK is smaller than the light cone sum rules estimation or the "tree" amplitude has been over estimated in QCDF. The latter possibility could happen if charming penguins are nonperturbative and not small, as argued in soft collinear effective theory. To differentiate between these two possibilities, we examine the similar B s → ρ + K − decay with the same technique. It is found that a large part of the uncertainties are canceled in the ratioIn QCDF, it is predicted to be 2.5 ± 0.2 which is independent on the form factor. However if charming penguins are important, this ratio could be very different from the QCDF prediction. Therefore the ratio of these two branching ratios could be an interesting indicator of the role of charming penguins in charmless B decays.