  We computed the leading order Wilson coefficients relevant to all the
exclusive $b\to s\ell^+\ell^-$ decays in the framework of Two Higgs Doublet
Model (2HDM) with a softly broken $\mathbb{Z}_2$ symmetry by including the
$\mathcal{O}(m_b)$ corrections. We elucidate the issue of appropriate matching
between the full and the effective theory when dealing with the (pseudo-)scalar
operators for which keeping the external momenta different from zero is
necessary. We then make a phenomenological analysis by using the measured
$\mathcal{B}(B_s\to \mu^+\mu^-)$ and $\mathcal{B}(B\to K
\mu^+\mu^-)_{\mathrm{high}-q^2}$, for which the hadronic uncertainties are well
controlled, and discussed their impact on various types of 2HDM. A brief
discussion of the decays with $\tau$-leptons in the final state is also made.
