Reynolds-averaged Navier-Stokes (RANS) equations are widely used in engineering turbulent flow simulations. However, RANS predictions may have large discrepancies due to the uncertainties in modeled Reynolds stresses. Recently, Wang et al. demonstrated that machine learning can be used to improve the RANS modeled Reynolds stresses by leveraging data from high fidelity simulations (Physics informed machine learning approach for reconstructing Reynolds stress modeling discrepancies based on DNS data. Physical Review Fluids. 2, 034603, 2017). However, solving for mean flows from the improved Reynolds stresses still poses significant challenges due to potential ill-conditioning of RANS equations with Reynolds stress closures. Enabling improved predictions of mean velocities are of profound practical importance, because often the velocity and its derived quantities (QoIs, e.g., drag, lift, surface friction), and not the Reynolds stress itself, are of ultimate interest in RANS simulations. To this end, we present a comprehensive framework for augmenting turbulence models with physics-informed machine learning, illustrating a complete workflow from identification of input features to final prediction of mean velocities. This work has two innovations. First, we demonstrate a systematic procedure to generate mean flow features based on the integrity basis for mean flow tensors. Second, we propose using machine learning to predict linear and nonlinear parts of the Reynolds stress tensor separately. Inspired by the finite polynomial representation of tensors in classical turbulence modeling, such a decomposition is instrumental in overcoming the ill-conditioning of RANS equations. Numerical tests demonstrated merits of the proposed framework.