In this paper we perform a comprehensive study of the four B → Kη (′) decays in the perturbative QCD (pQCD) factorization approach. We calculate the CP-averaged branching ratios and CP-violating asymmetries of B → Kη (′) decays in the ordinary η-η ′ , the η-η ′ -G, and the ηη ′ -G-η c mixing schemes. Besides the full leading order (LO) contributions, all currently known next-to-leading order (NLO) contributions to B → Kη (′) decays in the pQCD approach are taken into account. From our calculations and phenomenological analysis, we find the following. (a) The NLO contributions in general can provide significant enhancements to the LO pQCD predictions for the decay rates of the two B → Kη ′ decays, around 70%−89% in magnitude, but result in relatively small changes to Br(B → Kη). (b) Although the NLO pQCD predictions in all three considered mixing schemes agree well with the data within one standard deviation, those pQCD predictions in the η-η ′ -G mixing scheme provide the best interpretation for the measured pattern of Br(B → Kη (′) ): Br(B 0 → K 0 η) ≈ 1.13 × 10 −6 , Br(B 0 → K 0 η ′ ) ≈ 66.5 × 10 −6 , Br(B ± → K ± η) ≈ 2.36 × 10 −6 , and Br(B ± → K ± η ′ ) ≈ 67.3 × 10 −6 , which agree perfectly with the measured values. (c) The NLO pQCD predictions for the CP-violating asymmetries for the four considered decays are also in good consistent with the data. (d) The newly known NLO contribution to the relevant form factors M F F can produce about a 20% enhancement to the branching ratios Br(B → Kη ′ ), which plays an important role in closing the gap between the pQCD predictions and the relevant data. (e) The general expectations about the relative strength of the LO and NLO contributions from different sources are examined and confirmed by explicit numerical calculations.