We perform a global fit of the available polarized Semi-Inclusive Deep Inelastic Scattering (SIDIS), polarized pion-induced Drell-Yan (DY) and W-+/-/Z boson production data at (NLO)-L-3 and NNLO accuracy of the Transverse Momentum Dependent (TMD) evolution, and extract the Sivers function for u, d, s and for sea quarks. The Qiu-Sterman function is determined in a model independent way via the operator product expansion from the extracted Sivers function. The analysis is supplemented by additional studies, such as the estimation of applicability region, the impact of the unpolarized distributions' uncertainties, the universality of the Sivers functions, positivity constraints, the significance of the sign-change relation, and the comparison with the existing extractions.