We develop a numerical scheme for the calculation of tunneling current I and differential conductance dI/dV of metal- and CO-terminated STM tips on the topological insulators Bi2Se3, Bi2Te2Se as well as Bi2Te3, and find excellent agreement with experiment. The calculation is an application of Chen's derivative rule, whereby the Bloch functions are obtained from Wannier-interpolated tight-binding Hamiltonians and maximally localized Wannier functions from first-principle DFT+GW calculations. We observe signatures of the topological boundary modes, their hybridization with bulk bands, Van Hove singularities of the bulk bands, and characterize the orbital character of theses electronic modes using the high spatial resolution of STM and AFM. Bare DFT calculations are insufficient to explain the experimental data, which are instead accurately reproduced by many-body-corrected GW calculations.