We prove that the squared singular values of a fixed matrix multiplied with a truncation of a Haar distributed unitary matrix are distributed by a polynomial ensemble. This result is applied to a multiplication of a truncated unitary matrix with a random matrix. We show that the structure of polynomial ensembles and of certain Pfaffian ensembles is preserved. Furthermore we derive the joint singular value density of a product of truncated unitary matrices and its corresponding correlation kernel which can be written as a double contour integral. This leads to hard edge scaling limits that also include new finite rank perturbations of the Meijer G-kernels found for products of complex Ginibre random matrices.