In recent years, several experiments at the e−p collider HERA have collected high precision deepinelastic scattering (DIS) data on the spectrum of leading nucleon carrying a large fraction of the proton’s energy. In this paper, we have analyzed recent experimental data on the production of forward protons and neutrons in DIS at HERA in the framework of a perturbative QCD. We propose a technique based on the fractures functions framework, and extract the nucleon fracture functions (FFs) Mðn=pÞ 2 ðx;Q2; xLÞ from global QCD analysis of DIS data measured by the ZEUS Collaboration at HERA. We have shown that an approach based on the fracture functions formalism allows us to phenomenologically parametrize the nucleon FFs. Considering both leading neutron as well as leading proton production data at HERA, we present the results for the separate parton distributions for all parton species, including valence quark densities, the antiquark densities, the strange sea distribution, and the gluon distribution functions. We proposed several parametrizations for the nucleon FFs and open the possibility of these asymmetries. The obtained optimum set of nucleon FFs is accompanied by Hessian uncertainty sets which allow one to propagate uncertainties to other observables interest. The extracted results for the t-integrated leading neutron FLNð3Þ 2 ðx;Q2; xLÞ and leading proton FLPð3Þ 2 ðx;Q2; xLÞ structure functions are in good agreement with all data analyzed, for a wide range of fractional momentum variable x as well as the longitudinal momentum fraction xL.