The use of atomic orbitals in Hedin's GW approximation provides, in principle, an inexpensive alternative to plane-wave basis sets, especially when modeling large molecules. However, benchmarking of the algorithms and basis sets is essential for a careful balance between cost and accuracy. In this paper, we present an implementation of the GW approximation using numerical atomic orbitals and a pseudopotential treatment of core electrons. The combination of a contour deformation technique with a one-shot extraction of quasiparticle energies provides an efficient scheme for many applications. The performance of the implementation with respect to the basis set convergence and the effect of the use of pseudopotentials has been tested for the 117 closed-shell molecules from the G2/97 test set and 24 larger acceptor molecules from another recently proposed test set. Moreover, to demonstrate the potential of our method, we compute the thermally averaged GW density of states of a large photochromic compound by sampling ab initio molecular dynamics trajectories at different temperatures. The computed thermal line widths indicate approximately twice as large electron-phonon couplings with GW than with standard DFT-GGA calculations. This is further confirmed using frozen-phonon calculations.