We carry out a coarse-grained molecular dynamics simulation of phospholipid vesicles with transmembrane proteins. We measure the mean and Gaussian curvatures of our protein-embedded vesicles and quantitatively show how protein clusters change the shapes of their host vesicles. The effects of depletion force and vesiculation on protein clustering are also investigated. By increasing the protein concentration, clusters are fragmented to smaller bundles, which are then redistributed to form more symmetric structures corresponding to lower bending energies. Big clusters and highly aspherical vesicles cannot be formed when the fraction of protein to lipid molecules is large.