Computational fluid dynamics using the Reynolds-averaged Navier-Stokes (RANS) remains the most cost-effective approach to study wake flows and power losses in wind farms. The underlying assumptions associated with turbulence closures are the biggest sources of errors and uncertainties in the model predictions. This work aims to quantify model-form uncertainties in RANS simulations of wind farms at high Reynolds numbers under neutrally stratified conditions by perturbing the Reynolds stress tensor through a data-driven machine-learning technique. To this end, a two-step feature-selection method is applied to determine key features of the model. Then, the extreme gradient boosting algorithm is validated and employed to predict the perturbation amount and direction of the modeled Reynolds stress toward the limiting states of turbulence on the barycentric map. This procedure leads to a more accurate representation of the Reynolds stress anisotropy. The data-driven model is trained on high-fidelity data obtained from large-eddy simulation of a specific wind farm, and it is tested on two other (unseen) wind farms with distinct layouts to analyze its performance in cases with different turbine spacing and partial wake. The results indicate that, unlike the data-free approach in which a uniform and constant perturbation amount is applied to the entire computational domain, the proposed framework yields an optimal estimation of the uncertainty bounds for the RANS-predicted quantities of interest, including the wake velocity, turbulence intensity, and power losses in wind farms.