In this paper, we study a fully discrete vector penalty-projection method for the time-dependent incompressible magnetohydrodynamics flows. This fully discrete scheme is a combination of a mixed finite element approximation for spatial discretization and first-order backward Euler for temporal discretization. Moreover, unconditionally energy stable is established, and error estimates for the fully discrete scheme is also derived. Finally, some numerical experiments are provided to verify the theoretical results and illustrate the accuracy and efficiency of the proposed scheme.